Skip to content

P
pyMKcam
Commit 107ff4dc authored by sumpfralle's avatar sumpfralle
Browse files
Options

fixed code-style issues 


git-svn-id: https://pycam.svn.sourceforge.net/svnroot/pycam/trunk@491 bbaffbd6-741e-11dd-a85d-61de82d9cad9
parent 93ba4797
Hide whitespace changes
Inline Side-by-side
Showing with 742 additions and 646 deletions
src/pycam/Geometry/Line.py
View file @ 107ff4dc
......@@ -23,7 +23,7 @@ along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
try:
import OpenGL.GL as GL
GL_enabled = True
except:
except ImportError:
GL_enabled = False
......@@ -31,16 +31,20 @@ import math
class Line:
id=0
def __init__(self,p1,p2):
id = 0
def __init__(self, p1, p2):
self.id = Line.id
Line.id += 1
self.p1 = p1
self.p2 = p2
self._dir = None
self._len = None
def __repr__(self):
return "Line<%g,%g,%g>-<%g,%g,%g>" % (self.p1.x,self.p1.y,self.p1.z,
self.p2.x,self.p2.y,self.p2.z)
return "Line<%g,%g,%g>-<%g,%g,%g>" % (self.p1.x, self.p1.y, self.p1.z,
self.p2.x, self.p2.y, self.p2.z)
def __cmp__(self, other):
""" Two lines are equal if both pairs of points are at the same
locations.
......@@ -58,13 +62,13 @@ class Line:
return cmp(str(self), str(other))
def dir(self):
if not hasattr(self,"_dir"):
if self._dir is None:
self._dir = self.p2.sub(self.p1)
self._dir.normalize()
return self._dir
def len(self):
if not hasattr(self,"_len"):
if self._len is None:
self._len = self.p2.sub(self.p1).norm()
return self._len
......@@ -73,14 +77,14 @@ class Line:
def closest_point(self, p):
v = self.dir()
l = self.p1.dot(v)-p.dot(v)
l = self.p1.dot(v) - p.dot(v)
return self.p1.sub(v.mul(l))
def dist_to_point_sq(self, p):
return p.sub(self.closest_point(p)).normsq()
def dist_to_point(self, p):
return sqrt(self.dist_to_point_sq(p))
return math.sqrt(self.dist_to_point_sq(p))
def minx(self):
return min(self.p1.x, self.p2.x)
......@@ -117,10 +121,14 @@ class Line:
line_size = math.sqrt((line[0] ** 2) + (line[1] ** 2))
ortho_size = math.sqrt((ortho[0] ** 2) + (ortho[1] ** 2))
ortho_dest_size = line_size / 10.0
ortho = (ortho[0] * ortho_dest_size / ortho_size, ortho[1] * ortho_dest_size / ortho_size)
line_back = (-line[0] * ortho_dest_size / line_size, -line[1] * ortho_dest_size / line_size)
p3 = (self.p2.x + ortho[0] + line_back[0], self.p2.y + ortho[1] + line_back[1], self.p2.z)
p4 = (self.p2.x - ortho[0] + line_back[0], self.p2.y - ortho[1] + line_back[1], self.p2.z)
ortho = (ortho[0] * ortho_dest_size / ortho_size,
ortho[1] * ortho_dest_size / ortho_size)
line_back = (-line[0] * ortho_dest_size / line_size,
-line[1] * ortho_dest_size / line_size)
p3 = (self.p2.x + ortho[0] + line_back[0],
self.p2.y + ortho[1] + line_back[1], self.p2.z)
p4 = (self.p2.x - ortho[0] + line_back[0],
self.p2.y - ortho[1] + line_back[1], self.p2.z)
GL.glVertex3f(p3[0], p3[1], p3[2])
GL.glVertex3f(self.p2.x, self.p2.y, self.p2.z)
GL.glVertex3f(p4[0], p4[1], p4[2])
......
src/pycam/Geometry/Matrix.py
View file @ 107ff4dc
......@@ -20,13 +20,13 @@ You should have received a copy of the GNU General Public License
along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
"""
"""
various matrix related functions for PyCAM
"""
# various matrix related functions for PyCAM
from pycam.Geometry.Point import Point
import math
def get_dot_product(a, b):
""" calculate the dot product of two 3d vectors
......@@ -37,7 +37,7 @@ def get_dot_product(a, b):
@rtype: float
@return: the dot product is (a0*b0 + a1*b1 + a2*b2)
"""
return sum(map(lambda l1, l2: l1 * l2, a, b))
return sum(l1 * l2 for l1, l2 in zip(a, b))
def get_cross_product(a, b):
""" calculate the cross product of two 3d vectors
......@@ -103,7 +103,8 @@ def get_rotation_matrix_from_to(v_orig, v_dest):
arcsin = -1.0
rot_angle = math.asin(arcsin)
# calculate the rotation axis
# the rotation axis is equal to the cross product of the original and destination vectors
# The rotation axis is equal to the cross product of the original and
# destination vectors.
rot_axis = Point(v_orig[1] * v_dest[2] - v_orig[2] * v_dest[1],
v_orig[2] * v_dest[0] - v_orig[0] * v_dest[2],
v_orig[0] * v_dest[1] - v_orig[1] * v_dest[0])
......@@ -124,7 +125,7 @@ def get_rotation_matrix_from_to(v_orig, v_dest):
t * rot_axis.z * rot_axis.z + c)
def get_rotation_matrix_axis_angle(rot_axis, rot_angle):
""" calculate rotation matrix for a normalized "rot_axis" vector and an angle
""" calculate rotation matrix for a normalized vector and an angle
see http://mathworld.wolfram.com/RotationMatrix.html
@type rot_axis: tuple(float)
......
src/pycam/Geometry/Model.py
View file @ 107ff4dc
......@@ -25,13 +25,7 @@ import pycam.Exporters.STLExporter
from pycam.Geometry import Triangle, Line, Point
from pycam.Geometry.TriangleKdtree import TriangleKdtree
from pycam.Toolpath import Bounds
from utils import INFINITE
try:
import OpenGL.GL as GL
GL_enabled = True
except:
GL_enabled = False
from pycam.Geometry.utils import INFINITE
MODEL_TRANSFORMATIONS = {
......@@ -146,9 +140,12 @@ class BaseModel(object):
for point in item.get_points():
if not point.id in processed:
processed.append(point.id)
x = point.x * matrix[0][0] + point.y * matrix[0][1] + point.z * matrix[0][2] + matrix[0][3]
y = point.x * matrix[1][0] + point.y * matrix[1][1] + point.z * matrix[1][2] + matrix[1][3]
z = point.x * matrix[2][0] + point.y * matrix[2][1] + point.z * matrix[2][2] + matrix[2][3]
x = point.x * matrix[0][0] + point.y * matrix[0][1] \
+ point.z * matrix[0][2] + matrix[0][3]
y = point.x * matrix[1][0] + point.y * matrix[1][1] \
+ point.z * matrix[1][2] + matrix[1][3]
z = point.x * matrix[2][0] + point.y * matrix[2][1] \
+ point.z * matrix[2][2] + matrix[2][3]
point.x = x
point.y = y
point.z = z
......@@ -188,6 +185,7 @@ class Model(BaseModel):
self._kdtree_dirty = True
# enable/disable kdtree
self._use_kdtree = use_kdtree
self._t_kdtree = None
def append(self, item):
super(Model, self).append(item)
......@@ -313,7 +311,8 @@ class ContourModel(BaseModel):
finished = False
while not finished:
if len(new_group) > 1:
# calculate new intersections for each pair of adjacent lines
# Calculate new intersections for each pair of adjacent
# lines.
for index in range(len(new_group)):
if (index == 0) and (not closed_group):
# skip the first line if the group is not closed
......@@ -324,7 +323,8 @@ class ContourModel(BaseModel):
do_lines_intersection(l1, l2)
# Remove all lines that were marked as obsolete during
# intersection calculation.
clean_group = [line for line in new_group if not line.p1 is None]
clean_group = [line for line in new_group
if not line.p1 is None]
finished = len(new_group) == len(clean_group)
if (len(clean_group) == 1) and closed_group:
new_group = []
......
src/pycam/Geometry/Path.py
View file @ 107ff4dc
......@@ -26,6 +26,7 @@ don't really need complete "Point" instances that consume a lot of memory.
Since python 2.6 the "namedtuple" factory is available.
This reduces the memory consumption of a toolpath down to 1/3.
"""
try:
# this works for python 2.6 or above (saves memory)
import collections.namedtuple
......@@ -43,7 +44,7 @@ class Path:
Path.id += 1
self.top_join = None
self.bot_join = None
self.winding=0
self.winding = 0
self.points = []
def __repr__(self):
......@@ -54,7 +55,7 @@ class Path:
if first:
first = False
else:
s +="-"
s += "-"
s += "%d(%g,%g,%g)" % (p.id, p.x, p.y, p.z)
return s
......
src/pycam/Geometry/Plane.py
View file @ 107ff4dc
......@@ -20,7 +20,6 @@ You should have received a copy of the GNU General Public License
along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
"""
from Point import *
class Plane:
id = 0
......@@ -31,5 +30,5 @@ class Plane:
self.n = n
def __repr__(self):
return "Plane<%s,%s>" % (self.p,self.n)
return "Plane<%s,%s>" % (self.p, self.n)
src/pycam/Geometry/Point.py
View file @ 107ff4dc
......@@ -27,17 +27,19 @@ def _is_near(x, y):
class Point:
id=0
id = 0
def __init__(self,x,y,z):
def __init__(self, x, y, z):
self.id = Point.id
Point.id += 1
self.x = float(x)
self.y = float(y)
self.z = float(z)
self._norm = None
self._normsq = None
def __repr__(self):
return "Point%d<%g,%g,%g>" % (self.id,self.x,self.y,self.z)
return "Point%d<%g,%g,%g>" % (self.id, self.x, self.y, self.z)
def __cmp__(self, other):
""" Two points are equal if all dimensions are identical.
......@@ -57,30 +59,31 @@ class Point:
return cmp(str(self), str(other))
def mul(self, c):
return Point(self.x*c,self.y*c,self.z*c)
return Point(self.x * c, self.y * c, self.z * c)
def div(self, c):
return Point(self.x/c,self.y/c,self.z/c)
return Point(self.x / c, self.y / c, self.z / c)
def add(p1, p2):
return Point(p1.x+p2.x,p1.y+p2.y,p1.z+p2.z)
def add(self, p):
return Point(self.x + p.x, self.y + p.y, self.z + p.z)
def sub(p1, p2):
return Point(p1.x-p2.x,p1.y-p2.y,p1.z-p2.z)
def sub(self, p):
return Point(self.x - p.x, self.y - p.y, self.z - p.z)
def dot(p1, p2):
return p1.x*p2.x + p1.y*p2.y + p1.z*p2.z
def dot(self, p):
return self.x * p.x + self.y * p.y + self.z * p.z
def cross(p1, p2):
return Point(p1.y*p2.z-p2.y*p1.z, p2.x*p1.z-p1.x*p2.z, p1.x*p2.y-p2.x*p1.y)
def cross(self, p):
return Point(self.y * p.z - p.y * self.z, p.x * self.z - self.x * p.z,
self.x * p.y - p.x * self.y)
def normsq(self):
if not hasattr(self, "_normsq"):
if self._normsq is None:
self._normsq = self.dot(self)
return self._normsq
def norm(self):
if not hasattr(self, "_norm"):
if self._norm is None:
self._norm = math.sqrt(self.normsq())
return self._norm
......
src/pycam/Geometry/PointKdtree.py
View file @ 107ff4dc
......@@ -20,21 +20,21 @@ You should have received a copy of the GNU General Public License
along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
"""
import math
from Point import *
from Line import *
from Triangle import *
from kdtree import *
from pycam.Geometry.utils import epsilon
from pycam.Geometry.Point import Point
from pycam.Geometry.kdtree import Node, kdtree
class PointKdtree(kdtree):
def __init__(self, points=[],cutoff=5, cutoff_distance=0.5, tolerance=0.001):
def __init__(self, points=None, cutoff=5, cutoff_distance=0.5,
tolerance=epsilon):
if points is None:
points = []
self._n = None
self.tolerance=tolerance
self.tolerance = tolerance
nodes = []
for p in points:
n = Node();
n = Node()
n.point = p
n.bound = []
n.bound.append(p.x)
......@@ -54,19 +54,18 @@ class PointKdtree(kdtree):
if self._n:
n = self._n
else:
n = Node();
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:
(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)
n.p = Point(x, y, z)
self._n = None
self.insert(n)
return n.p
src/pycam/Geometry/PolygonExtractor.py
View file @ 107ff4dc
......@@ -20,21 +20,21 @@ You should have received a copy of the GNU General Public License
along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
"""
from pycam.Utils.iterators import *
from pycam.Utils.iterators import Iterator, CyclicIterator
from pycam.Geometry.utils import *
from pycam.Geometry.Path import *
from pycam.Geometry.Point import *
from pycam.Geometry.utils import INFINITE
from pycam.Geometry.Path import Path
from pycam.Exporters.SVGExporter import SVGExporter
DEBUG_POLYGONEXTRACTOR=False
DEBUG_POLYGONEXTRACTOR2=False
DEBUG_POLYGONEXTRACTOR3=False
DEBUG_POLYGONEXTRACTOR = False
DEBUG_POLYGONEXTRACTOR2 = False
DEBUG_POLYGONEXTRACTOR3 = False
class PolygonExtractor:
CONTOUR=1
MONOTONE=2
CONTOUR = 1
MONOTONE = 2
def __init__(self, policy=MONOTONE):
self.policy = policy
......@@ -42,20 +42,21 @@ class PolygonExtractor:
self.ver_path_list = None
self.merge_path_list = None
def new_direction(self, dir):
self.current_dir = dir
def new_direction(self, direction):
self.current_dir = direction
self.all_path_list = []
self.curr_path_list = []
self.prev_line = []
self.curr_line = []
if DEBUG_POLYGONEXTRACTOR3:
self.svg = SVGExporter("test-%d.svg" % dir)
if dir == 0:
self.svg = SVGExporter("test-%d.svg" % direction)
if direction == 0:
self.svg.fill("red")
else:
self.svg.fill("blue")
if self.policy == PolygonExtractor.CONTOUR and dir == 1 and self.hor_path_list:
if (self.policy == PolygonExtractor.CONTOUR) and (direction == 1) \
and self.hor_path_list:
self.last_x = -INFINITE
self.delta_x = 0
self.convert_hor_path_list()
......@@ -75,63 +76,71 @@ class PolygonExtractor:
if self.policy == PolygonExtractor.CONTOUR and self.hor_path_list:
for path in self.all_path_list:
if DEBUG_POLYGONEXTRACTOR2: print "points=", path.points
if DEBUG_POLYGONEXTRACTOR2:
print "points=", path.points
i = 0
while i < len(path.points)-1:
if path.points[i].x > path.points[i+1].x:
if DEBUG_POLYGONEXTRACTOR2: print "drop point %d:" % path.points[i].id
if DEBUG_POLYGONEXTRACTOR2:
print "drop point %d:" % path.points[i].id
path.points = path.points[:i] + path.points[i+1:]
if i>0:
if i > 0:
i -= 1
else:
i += 1
if DEBUG_POLYGONEXTRACTOR: print "%d paths" % len(self.all_path_list)
for path in self.all_path_list:
if DEBUG_POLYGONEXTRACTOR: print "%d:" % path.id,
if DEBUG_POLYGONEXTRACTOR: print "%d ->" % path.top_join.id
for point in path.points:
if DEBUG_POLYGONEXTRACTOR: print "%d(%g,%g)" % (point.id, point.x, point.y),
if DEBUG_POLYGONEXTRACTOR: print "->%d" % path.bot_join.id
if DEBUG_POLYGONEXTRACTOR:
print "%d paths" % len(self.all_path_list)
for path in self.all_path_list:
print "%d:" % path.id,
print "%d ->" % path.top_join.id
for point in path.points:
print "%d(%g,%g)" % (point.id, point.x, point.y),
print "->%d" % path.bot_join.id
path_list = []
while len(self.all_path_list)>0:
while len(self.all_path_list) > 0:
p0 = self.all_path_list[0]
for path in self.all_path_list:
if path.id < p0.id:
p0 = path
if DEBUG_POLYGONEXTRACTOR: print "linking %d" % p0.id
if DEBUG_POLYGONEXTRACTOR:
print "linking %d" % p0.id
self.all_path_list.remove(p0)
p1 = p0.bot_join
while True:
if DEBUG_POLYGONEXTRACTOR: print "splice %d into %d" % (p1.id, p0.id)
if DEBUG_POLYGONEXTRACTOR:
print "splice %d into %d" % (p1.id, p0.id)
self.all_path_list.remove(p1)
p1.reverse()
p0.points += p1.points
if p1.top_join == p0:
break;
break
p2 = p1.top_join
if DEBUG_POLYGONEXTRACTOR: print "splicing %d into %d" % (p2.id, p0.id)
if DEBUG_POLYGONEXTRACTOR:
print "splicing %d into %d" % (p2.id, p0.id)
self.all_path_list.remove(p2)
p0.points += p2.points
p1 = p2.bot_join
path_list.append(p0)
if DEBUG_POLYGONEXTRACTOR: print "%d paths" % len(path_list)
for path in path_list:
if DEBUG_POLYGONEXTRACTOR: print "path %d(w=%d): " % (path.id, path.winding),
for point in path.points:
if DEBUG_POLYGONEXTRACTOR: print "%d(%g,%g)" % (point.id, point.x, point.y),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR:
print "%d paths" % len(path_list)
for path in path_list:
print "path %d(w=%d): " % (path.id, path.winding),
for point in path.points:
print "%d(%g,%g)" % (point.id, point.x, point.y),
print
if self.current_dir==0:
if self.current_dir == 0:
self.hor_path_list = path_list
elif self.current_dir==1:
if self.policy == PolygonExtractor.CONTOUR and self.hor_path_list and path_list:
elif self.current_dir == 1:
if (self.policy == PolygonExtractor.CONTOUR) \
and self.hor_path_list and path_list:
self.merge_path_list = path_list
else:
self.ver_path_list = path_list
......@@ -172,12 +181,12 @@ class PolygonExtractor:
self.curr_line.append(p)
def end_scanline(self):
if self.current_dir==0:
if self.current_dir == 0:
self.process_hor_scanline(self.curr_line)
elif self.current_dir==1:
elif self.current_dir == 1:
if self.policy == PolygonExtractor.CONTOUR and self.hor_path_list:
next_x = -INFINITE
if len(self.curr_line)>0:
if len(self.curr_line) > 0:
next_x = self.curr_line[0].x
self.delta_x = next_x - self.last_x
self.last_x = next_x
......@@ -185,99 +194,112 @@ class PolygonExtractor:
next_x = self.last_x + self.delta_x
if next_x > -INFINITE:
self.process_virtual_hor_scanline(next_x)
if DEBUG_POLYGONEXTRACTOR2:
if DEBUG_POLYGONEXTRACTOR2:
print "scanline =", self.curr_line
self.process_ver_scanline(self.curr_line)
def process_hor_scanline(self, scanline):
if DEBUG_POLYGONEXTRACTOR:
if DEBUG_POLYGONEXTRACTOR:
last = 0
inside = False
s = ""
for point in scanline:
next = point.x
next_x = point.x
if inside:
s += "*" * int(next-last)
s += "*" * int(next_x - last)
else:
s += " " * int(next-last)
last = next
s += " " * int(next_x - last)
last = next_x
inside = not inside
print s
if DEBUG_POLYGONEXTRACTOR: print "active paths: ",
for path in self.curr_path_list:
if DEBUG_POLYGONEXTRACTOR: print "%d(%g,%g)" % (path.id, path.points[-1].x, path.points[-1].y),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR: print "prev points: ",
for point in self.prev_line:
if DEBUG_POLYGONEXTRACTOR: print "(%g,%g)" % (point.x, point.y),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR: print "active points: ",
for point in scanline:
if DEBUG_POLYGONEXTRACTOR: print "%d(%g,%g)" % (point.id, point.x, point.y),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR:
print "active paths: ",
for path in self.curr_path_list:
print "%d(%g,%g)" \
% (path.id, path.points[-1].x, path.points[-1].y),
print
print "prev points: ",
for point in self.prev_line:
print "(%g,%g)" % (point.x, point.y),
print
print "active points: ",
for point in scanline:
print "%d(%g,%g)" % (point.id, point.x, point.y),
print
prev_point = Iterator(self.prev_line)
curr_point = Iterator(scanline)
curr_path = Iterator(self.curr_path_list)
winding = 0
while prev_point.remains()>0 or curr_point.remains()>0:
if DEBUG_POLYGONEXTRACTOR: print "num_prev=", prev_point.remains(), ", num_curr=",curr_point.remains()
if prev_point.remains()==0 and curr_point.remains()>=2:
while (prev_point.remains() > 0) or (curr_point.remains() > 0):
if DEBUG_POLYGONEXTRACTOR:
print "num_prev=%d, num_curr=%d" \
% (prev_point.remains(), curr_point.remains())
if (prev_point.remains() == 0) and (curr_point.remains() >= 2):
c0 = curr_point.next()
c1 = curr_point.next()
# new path starts
p0 = Path()
p0.winding = winding+1
if DEBUG_POLYGONEXTRACTOR: print "new path %d(%g,%g)" % ( p0.id, c0.x, c0.y)
p0.winding = winding + 1
if DEBUG_POLYGONEXTRACTOR:
print "new path %d(%g,%g)" % (p0.id, c0.x, c0.y)
p0.append(c0)
self.curr_path_list.append(p0)
p1 = Path()
p1.winding = winding
if DEBUG_POLYGONEXTRACTOR: print "new path %d(%g,%g)" % (p1.id, c1.x, c1.y)
if DEBUG_POLYGONEXTRACTOR:
print "new path %d(%g,%g)" % (p1.id, c1.x, c1.y)
p1.append(c1)
self.curr_path_list.append(p1)
p0.top_join = p1
p1.top_join = p0
continue
if prev_point.remains()>=2 and curr_point.remains()==0:
if (prev_point.remains() >= 2) and (curr_point.remains() == 0):
#old path ends
p0 = curr_path.takeNext()
if DEBUG_POLYGONEXTRACTOR: print "end path %d" % p0.id
if DEBUG_POLYGONEXTRACTOR:
print "end path %d" % p0.id
self.all_path_list.append(p0)
prev_point.next()
p1 = curr_path.takeNext()
if DEBUG_POLYGONEXTRACTOR: print "end path %d" % p1.id
if DEBUG_POLYGONEXTRACTOR:
print "end path %d" % p1.id
self.all_path_list.append(p1)
prev_point.next()
p0.bot_join = p1
p1.bot_join = p0
continue
if prev_point.remains()>=2 and curr_point.remains()>=2:
if (prev_point.remains() >= 2) and (curr_point.remains() >= 2):
p0 = prev_point.peek(0)
p1 = prev_point.peek(1)
c0 = curr_point.peek(0)
c1 = curr_point.peek(1)
if DEBUG_POLYGONEXTRACTOR: print "overlap test: p0=%g p1=%g" % (p0.x, p1.x)
if DEBUG_POLYGONEXTRACTOR: print "overlap test: c0=%g c1=%g" % (c0.x, c1.x)
if DEBUG_POLYGONEXTRACTOR:
print "overlap test: p0=%g p1=%g" % (p0.x, p1.x)
print "overlap test: c0=%g c1=%g" % (c0.x, c1.x)
if c1.x < p0.x:
# new segment is completely to the left
# new path starts
s0 = Path()
if DEBUG_POLYGONEXTRACTOR: print "new path %d(%g,%g) w=%d" % (s0.id, c0.x, c0.y, winding+1)
if DEBUG_POLYGONEXTRACTOR:
print "new path %d(%g,%g) w=%d" \
% (s0.id, c0.x, c0.y, winding + 1)
s0.append(c0)
curr_path.insert(s0)
s1 = Path()
s0.winding = winding+1
s0.winding = winding + 1
s1.winding = winding
if DEBUG_POLYGONEXTRACTOR: print "new path %d(%g,%g) w=%d" % (s1.id, c1.x, c1.y, winding)
if DEBUG_POLYGONEXTRACTOR:
print "new path %d(%g,%g) w=%d" \
% (s1.id, c1.x, c1.y, winding)
s1.append(c1)
curr_path.insert(s1)
curr_point.next()
......@@ -288,16 +310,18 @@ class PolygonExtractor:
# new segment is completely to the right
# old path ends
s0 = curr_path.takeNext()
if DEBUG_POLYGONEXTRACTOR: print "end path %d" % s0.id
if DEBUG_POLYGONEXTRACTOR:
print "end path %d" % s0.id
self.all_path_list.append(s0)
prev_point.next()
s1 = curr_path.takeNext()
if DEBUG_POLYGONEXTRACTOR: print "end path %d" % s1.id
if DEBUG_POLYGONEXTRACTOR:
print "end path %d" % s1.id
self.all_path_list.append(s1)
prev_point.next()
s0.bot_join = s1;
s1.bot_join = s0;
winding = s1.winding;
s0.bot_join = s1
s1.bot_join = s0
winding = s1.winding
else:
# new segment is overlapping
left_path = curr_path.next()
......@@ -316,8 +340,10 @@ class PolygonExtractor:
# check for path joins
if prev_point.remains()>=2:
p2 = prev_point.peek(1)
if DEBUG_POLYGONEXTRACTOR: print "join test: p0=%g p1=%g p2=%g" % (p0.x, p1.x, p2.x)
if DEBUG_POLYGONEXTRACTOR: print "join test: c0=%g c1=%g" % (c0.x, c1.x)
if DEBUG_POLYGONEXTRACTOR:
print "join test: p0=%g p1=%g p2=%g" \
% (p0.x, p1.x, p2.x)
print "join test: c0=%g c1=%g" % (c0.x, c1.x)
if p2.x <= c1.x:
overlap_p = True
if self.policy == PolygonExtractor.CONTOUR:
......@@ -326,15 +352,19 @@ class PolygonExtractor:
if curr_path.remains()>=1:
right_path = curr_path.peek()
self.all_path_list.append(s0)
self.all_path_list.append(s1);
if DEBUG_POLYGONEXTRACTOR: print "path %d joins %d" % (s0.id, s1.id)
s0.bot_join = s1;
s1.bot_join = s0;
self.all_path_list.append(s1)
if DEBUG_POLYGONEXTRACTOR:
print "path %d joins %d" \
% (s0.id, s1.id)
s0.bot_join = s1
s1.bot_join = s0
elif self.policy == PolygonExtractor.MONOTONE:
s0 = curr_path.takeNext()
left_path.bot_join = s0
s0.bot_join = left_path
if DEBUG_POLYGONEXTRACTOR: print "path %d joins %d" % (left_path.id, s0.id)
if DEBUG_POLYGONEXTRACTOR:
print "path %d joins %d" \
% (left_path.id, s0.id)
curr_path.remove(left_path)
self.all_path_list.append(left_path)
self.all_path_list.append(s0)
......@@ -353,17 +383,21 @@ class PolygonExtractor:
# check for path splits
if curr_point.remains()>=2:
c2 = curr_point.peek(1)
if DEBUG_POLYGONEXTRACTOR: print "split test: p0=%g p1=%g" % (p0.x, p1.x)
if DEBUG_POLYGONEXTRACTOR: print "split test: c0=%g c1=%g c2=%g" % (c0.x, c1.x, c2.x)
if DEBUG_POLYGONEXTRACTOR:
print "split test: p0=%g p1=%g" % (p0.x, p1.x)
print "split test: c0=%g c1=%g c2=%g" \
% (c0.x, c1.x, c2.x)
if c2.x <= p1.x:
overlap_c = True
s0 = Path()
s1 = Path()
s0.winding=winding+1
s1.winding=winding
s0.winding = winding + 1
s1.winding = winding
s0.top_join = s1
s1.top_join = s0
if DEBUG_POLYGONEXTRACTOR: print "region split into %d and %d (w=%d)" %(s0.id, s1.id, winding+1)
if DEBUG_POLYGONEXTRACTOR:
print "region split into %d and %d (w=%d)" \
% (s0.id, s1.id, winding + 1)
curr_point.next()
c0 = curr_point.next()
if self.policy == PolygonExtractor.CONTOUR:
......@@ -377,26 +411,32 @@ class PolygonExtractor:
curr_path.insertBefore(s0)
curr_path.insertBefore(s1)
left_point = c2
if curr_point.remains()>=1:
if curr_point.remains() >= 1:
c1 = curr_point.peek(0)
right_point = c1
else:
overlap_c = False
if DEBUG_POLYGONEXTRACTOR: print "add to path %d(%g,%g)" % (left_path.id, left_point.x, left_point.y)
if DEBUG_POLYGONEXTRACTOR:
print "add to path %d(%g,%g)" \
% (left_path.id, left_point.x, left_point.y)
left_path.append(left_point)
right_path.append(right_point)
if right_path == curr_path.peek():
curr_path.next()
if DEBUG_POLYGONEXTRACTOR: print "add to path %d(%g,%g)" % (right_path.id, right_point.x, right_point.y)
winding = right_path.winding;
if DEBUG_POLYGONEXTRACTOR:
print "add to path %d(%g,%g)" \
% (right_path.id, right_point.x, right_point.y)
winding = right_path.winding
prev_point.next()
curr_point.next()
if DEBUG_POLYGONEXTRACTOR: print "active paths: ",
for path in self.curr_path_list:
if DEBUG_POLYGONEXTRACTOR: print "%d(%g,%g,w=%d)" % (path.id, path.points[-1].x, path.points[-1].y, path.winding),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR:
print "active paths: ",
for path in self.curr_path_list:
print "%d(%g,%g,w=%d)" % (path.id, path.points[-1].x,
path.points[-1].y, path.winding),
print
self.prev_line = scanline
......@@ -408,102 +448,114 @@ class PolygonExtractor:
self.cont.fill("red")
else:
self.cont.fill("blue")
self.cont.AddDot(p.x,p.y)
self.cont.AddDot(p.x, p.y)
self.cont.fill("black")
self.cont.AddText(p.x,p.y, str(p.id))
self.cont.AddText(p.x, p.y, str(p.id))
if prev:
self.cont.AddLine(prev.x,prev.y,p.x,p.y)
self.cont.AddLine(prev.x, prev.y, p.x, p.y)
prev = p
if DEBUG_POLYGONEXTRACTOR:
last = 0
inside = False
s = ""
for point in scanline:
next = point.y
next_y = point.y
if inside:
s += "*" * int(next-last)
s += "*" * int(next_y - last)
else:
s += " " * int(next-last)
last = next
s += " " * int(next_y - last)
last = next_y
inside = not inside
print s
if DEBUG_POLYGONEXTRACTOR: print "active paths: ",
for path in self.curr_path_list:
if DEBUG_POLYGONEXTRACTOR: print "%d(%g,%g)" % (path.id, path.points[-1].x, path.points[-1].y),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR: print "prev points: ",
for point in self.prev_line:
if DEBUG_POLYGONEXTRACTOR: print "(%g,%g)" % (point.x, point.y),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR: print "active points: ",
for point in scanline:
if DEBUG_POLYGONEXTRACTOR: print "%d(%g,%g)" % (point.id, point.x, point.y),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR:
print "active paths: ",
for path in self.curr_path_list:
print "%d(%g,%g)" \
% (path.id, path.points[-1].x, path.points[-1].y),
print
print "prev points: ",
for point in self.prev_line:
print "(%g,%g)" % (point.x, point.y),
print
print "active points: ",
for point in scanline:
print "%d(%g,%g)" % (point.id, point.x, point.y),
print
prev_point = Iterator(self.prev_line)
curr_point = Iterator(scanline)
curr_path = Iterator(self.curr_path_list)
winding = 0
while prev_point.remains()>0 or curr_point.remains()>0:
if DEBUG_POLYGONEXTRACTOR: print "num_prev=", prev_point.remains(), ", num_curr=",curr_point.remains()
if prev_point.remains()==0 and curr_point.remains()>=2:
while (prev_point.remains() > 0) or (curr_point.remains() > 0):
if DEBUG_POLYGONEXTRACTOR:
print "num_prev=%d, num_curr=%d" \
% (prev_point.remains(), curr_point.remains())
if (prev_point.remains() == 0) and (curr_point.remains() >= 2):
c0 = curr_point.next()
c1 = curr_point.next()
# new path starts
p0 = Path()
p0.winding = winding+1
if DEBUG_POLYGONEXTRACTOR: print "new path %d(%g,%g)" % ( p0.id, c0.x, c0.y)
p0.winding = winding + 1
if DEBUG_POLYGONEXTRACTOR:
print "new path %d(%g,%g)" % (p0.id, c0.x, c0.y)
p0.append(c0)
self.curr_path_list.append(p0)
p1 = Path()
p1.winding = winding
if DEBUG_POLYGONEXTRACTOR: print "new path %d(%g,%g)" % (p1.id, c1.x, c1.y)
if DEBUG_POLYGONEXTRACTOR:
print "new path %d(%g,%g)" % (p1.id, c1.x, c1.y)
p1.append(c1)
self.curr_path_list.append(p1)
p0.top_join = p1
p1.top_join = p0
continue
if prev_point.remains()>=2 and curr_point.remains()==0:
if (prev_point.remains() >= 2) and (curr_point.remains() == 0):
#old path ends
p0 = curr_path.takeNext()
if DEBUG_POLYGONEXTRACTOR: print "end path %d" % p0.id
if DEBUG_POLYGONEXTRACTOR:
print "end path %d" % p0.id
self.all_path_list.append(p0)
prev_point.next()
p1 = curr_path.takeNext()
if DEBUG_POLYGONEXTRACTOR: print "end path %d" % p1.id
if DEBUG_POLYGONEXTRACTOR:
print "end path %d" % p1.id
self.all_path_list.append(p1)
prev_point.next()
p0.bot_join = p1
p1.bot_join = p0
continue
if prev_point.remains()>=2 and curr_point.remains()>=2:
if (prev_point.remains() >= 2) and (curr_point.remains() >= 2):
p0 = prev_point.peek(0)
p1 = prev_point.peek(1)
c0 = curr_point.peek(0)
c1 = curr_point.peek(1)
if DEBUG_POLYGONEXTRACTOR: print "overlap test: p0=%g p1=%g" % (p0.x, p1.x)
if DEBUG_POLYGONEXTRACTOR: print "overlap test: c0=%g c1=%g" % (c0.x, c1.x)
if DEBUG_POLYGONEXTRACTOR:
print "overlap test: p0=%g p1=%g" % (p0.x, p1.x)
print "overlap test: c0=%g c1=%g" % (c0.x, c1.x)
if c1.y < p0.y:
# new segment is completely to the left
# new path starts
s0 = Path()
if DEBUG_POLYGONEXTRACTOR: print "new path %d(%g,%g) w=%d" % (s0.id, c0.x, c0.y, winding+1)
if DEBUG_POLYGONEXTRACTOR:
print "new path %d(%g,%g) w=%d" \
% (s0.id, c0.x, c0.y, winding + 1)
s0.append(c0)
curr_path.insert(s0)
s1 = Path()
s0.winding = winding+1
s0.winding = winding + 1
s1.winding = winding
if DEBUG_POLYGONEXTRACTOR: print "new path %d(%g,%g) w=%d" % (s1.id, c1.x, c1.y, winding)
if DEBUG_POLYGONEXTRACTOR:
print "new path %d(%g,%g) w=%d" \
% (s1.id, c1.x, c1.y, winding)
s1.append(c1)
curr_path.insert(s1)
curr_point.next()
......@@ -514,16 +566,18 @@ class PolygonExtractor:
# new segment is completely to the right
# old path ends
s0 = curr_path.takeNext()
if DEBUG_POLYGONEXTRACTOR: print "end path %d" % s0.id
if DEBUG_POLYGONEXTRACTOR:
print "end path %d" % s0.id
self.all_path_list.append(s0)
prev_point.next()
s1 = curr_path.takeNext()
if DEBUG_POLYGONEXTRACTOR: print "end path %d" % s1.id
if DEBUG_POLYGONEXTRACTOR:
print "end path %d" % s1.id
self.all_path_list.append(s1)
prev_point.next()
s0.bot_join = s1;
s1.bot_join = s0;
winding = s1.winding;
s0.bot_join = s1
s1.bot_join = s0
winding = s1.winding
else:
# new segment is overlapping
left_path = curr_path.next()
......@@ -540,27 +594,33 @@ class PolygonExtractor:
overlap_p = False
overlap_c = False
# check for path joins
if prev_point.remains()>=2:
if prev_point.remains() >= 2:
p2 = prev_point.peek(1)
if DEBUG_POLYGONEXTRACTOR: print "join test: p0=%g p1=%g p2=%g" % (p0.x, p1.x, p2.x)
if DEBUG_POLYGONEXTRACTOR: print "join test: c0=%g c1=%g" % (c0.x, c1.x)
if DEBUG_POLYGONEXTRACTOR:
print "join test: p0=%g p1=%g p2=%g" \
% (p0.x, p1.x, p2.x)
print "join test: c0=%g c1=%g" % (c0.x, c1.x)
if p2.y <= c1.y:
overlap_p = True
if self.policy == PolygonExtractor.CONTOUR:
s0 = curr_path.takeNext()
s1 = curr_path.takeNext()
if curr_path.remains()>=1:
if curr_path.remains() >= 1:
right_path = curr_path.peek()
self.all_path_list.append(s0)
self.all_path_list.append(s1);
if DEBUG_POLYGONEXTRACTOR: print "path %d joins %d" % (s0.id, s1.id)
s0.bot_join = s1;
s1.bot_join = s0;
self.all_path_list.append(s1)
if DEBUG_POLYGONEXTRACTOR:
print "path %d joins %d" \
% (s0.id, s1.id)
s0.bot_join = s1
s1.bot_join = s0
elif self.policy == PolygonExtractor.MONOTONE:
s0 = curr_path.takeNext()
left_path.bot_join = s0
s0.bot_join = left_path
if DEBUG_POLYGONEXTRACTOR: print "path %d joins %d" % (left_path.id, s0.id)
if DEBUG_POLYGONEXTRACTOR:
print "path %d joins %d" \
% (left_path.id, s0.id)
curr_path.remove(left_path)
self.all_path_list.append(left_path)
self.all_path_list.append(s0)
......@@ -571,7 +631,7 @@ class PolygonExtractor:
prev_point.next()
winding = s1.winding
p0 = p2
if prev_point.remains()>=1:
if prev_point.remains() >= 1:
p1 = prev_point.peek(0)
else:
overlap_p = False
......@@ -579,17 +639,21 @@ class PolygonExtractor:
# check for path splits
if curr_point.remains()>=2:
c2 = curr_point.peek(1)
if DEBUG_POLYGONEXTRACTOR: print "split test: p0=%g p1=%g" % (p0.x, p1.x)
if DEBUG_POLYGONEXTRACTOR: print "split test: c0=%g c1=%g c2=%g" % (c0.x, c1.x, c2.x)
if DEBUG_POLYGONEXTRACTOR:
print "split test: p0=%g p1=%g" % (p0.x, p1.x)
print "split test: c0=%g c1=%g c2=%g" \
% (c0.x, c1.x, c2.x)
if c2.y <= p1.y:
overlap_c = True
s0 = Path()
s1 = Path()
s0.winding=winding+1
s1.winding=winding
s0.winding = winding + 1
s1.winding = winding
s0.top_join = s1
s1.top_join = s0
if DEBUG_POLYGONEXTRACTOR: print "region split into %d and %d (w=%d)" %(s0.id, s1.id, winding+1)
if DEBUG_POLYGONEXTRACTOR:
print "region split into %d and %d (w=%d)" \
% (s0.id, s1.id, winding + 1)
curr_point.next()
c0 = curr_point.next()
if self.policy == PolygonExtractor.CONTOUR:
......@@ -603,31 +667,38 @@ class PolygonExtractor:
curr_path.insertBefore(s0)
curr_path.insertBefore(s1)
left_point = c2
if curr_point.remains()>=1:
if curr_point.remains() >= 1:
c1 = curr_point.peek(0)
right_point = c1
else:
overlap_c = False
if DEBUG_POLYGONEXTRACTOR: print "add to path %d(%g,%g)" % (left_path.id, left_point.x, left_point.y)
if DEBUG_POLYGONEXTRACTOR:
print "add to path %d(%g,%g)" \
% (left_path.id, left_point.x, left_point.y)
left_path.append(left_point)
right_path.append(right_point)
if right_path == curr_path.peek():
curr_path.next()
if DEBUG_POLYGONEXTRACTOR: print "add to path %d(%g,%g)" % (right_path.id, right_point.x, right_point.y)
winding = right_path.winding;
if DEBUG_POLYGONEXTRACTOR:
print "add to path %d(%g,%g)" \
% (right_path.id, right_point.x, right_point.y)
winding = right_path.winding
prev_point.next()
curr_point.next()
if DEBUG_POLYGONEXTRACTOR: print "active paths: ",
for path in self.curr_path_list:
if DEBUG_POLYGONEXTRACTOR: print "%d(%g,%g,w=%d)" % (path.id, path.points[-1].x, path.points[-1].y, path.winding),
if DEBUG_POLYGONEXTRACTOR: print
if DEBUG_POLYGONEXTRACTOR:
print "active paths: ",
for path in self.curr_path_list:
print "%d(%g,%g,w=%d)" % (path.id, path.points[-1].x,
path.points[-1].y, path.winding),
print
self.prev_line = scanline
def convert_hor_path_list(self):
if DEBUG_POLYGONEXTRACTOR2: print "converting hor paths"
if DEBUG_POLYGONEXTRACTOR2:
print "converting hor paths"
hor_path_list = []
for s in self.hor_path_list:
allsame = True
......@@ -635,14 +706,14 @@ class PolygonExtractor:
maxy = s.points[0].y
for p in s.points:
if not p.x == s.points[0].x:
allsame=False
allsame = False
if p.y < miny:
miny = p.y
if p.y > maxy:
maxy = p.y
if allsame:
if DEBUG_POLYGONEXTRACTOR2: print "all same !"
if DEBUG_POLYGONEXTRACTOR2:
print "all same !"
s0 = Path()
for p in s.points:
if p.y == miny:
......@@ -655,48 +726,52 @@ class PolygonExtractor:
hor_path_list.append(s1)
continue
prev = s.points[-1]
iter = CyclicIterator(s.points)
p_iter = CyclicIterator(s.points)
p = s.points[0]
next = iter.next()
while not (prev.x >= p.x and next.x > p.x):
p = next
next = iter.next()
next_p = p_iter.next()
while not ((prev.x >= p.x) and (next_p.x > p.x)):
p = next_p
next_p = p_iter.next()
count = 0
while count < len(s.points):
s0 = Path()
while next.x >= p.x:
while next_p.x >= p.x:
s0.append(p)
p = next
next = iter.next()
p = next_p
next_p = p_iter.next()
count += 1
s0.append(p)
while len(s0.points)>1 and s0.points[0].x == s0.points[1].x:
while (len(s0.points) > 1) \
and (s0.points[0].x == s0.points[1].x):
s0.points = s0.points[1:]
while len(s0.points)>1 and s0.points[-2].x == s0.points[-1].x:
while (len(s0.points) > 1) \
and (s0.points[-2].x == s0.points[-1].x):
s0.points = s0.points[0:-1]
hor_path_list.append(s0)
s1 = Path()
while next.x <= p.x:
while next_p.x <= p.x:
s1.append(p)
p = next
next = iter.next()
p = next_p
next_p = p_iter.next()
count += 1
s1.append(p)
s1.reverse()
while len(s1.points)>1 and s1.points[0].x == s1.points[1].x:
while (len(s1.points) > 1) \
and (s1.points[0].x == s1.points[1].x):
s1.points = s1.points[1:]
while len(s1.points)>1 and s1.points[-2].x == s1.points[-1].x:
s1.points = s1.points[0:-1]
while (len(s1.points) > 1) \
and (s1.points[-2].x == s1.points[-1].x):
s1.points = s1.points[:-1]
hor_path_list.append(s1)
hor_path_list.sort(cmp=lambda a,b:cmp(a.points[0].x,b.points[0].x))
if DEBUG_POLYGONEXTRACTOR2:
hor_path_list.sort(cmp=lambda a, b: cmp(a.points[0].x, b.points[0].x))
if DEBUG_POLYGONEXTRACTOR2:
print "ver_hor_path_list = ", hor_path_list
for s in hor_path_list:
print "s%d =" % s.id,
for point in s.points:
print point.id,
print
print
self.ver_hor_path_list = hor_path_list
self.act_hor_path_list = []
......@@ -709,50 +784,64 @@ class PolygonExtractor:
while next_x <= _next_x:
next_x = INFINITE
if len(self.ver_hor_path_list)>0 and self.ver_hor_path_list[0].points[0].x < next_x:
if self.ver_hor_path_list \
and (self.ver_hor_path_list[0].points[0].x < next_x):
next_x = self.ver_hor_path_list[0].points[0].x
if len(self.act_hor_path_list)>0 and self.act_hor_path_list[0].points[0].x < next_x:
if self.act_hor_path_list \
and (self.act_hor_path_list[0].points[0].x < next_x):
next_x = self.act_hor_path_list[0].points[0].x
if next_x >= _next_x:
return
if DEBUG_POLYGONEXTRACTOR2:
if DEBUG_POLYGONEXTRACTOR2:
print "ver_hor_path_list =", self.ver_hor_path_list
print "act_hor_path_list =", self.act_hor_path_list
print "next_x =", next_x
if len(self.ver_hor_path_list)>0 and self.ver_hor_path_list[0].points[0].x <= next_x:
while len(self.ver_hor_path_list)>0 and self.ver_hor_path_list[0].points[0].x <= next_x:
if self.ver_hor_path_list \
and (self.ver_hor_path_list[0].points[0].x <= next_x):
while self.ver_hor_path_list \
and (self.ver_hor_path_list[0].points[0].x <= next_x):
self.act_hor_path_list.append(self.ver_hor_path_list[0])
self.ver_hor_path_list = self.ver_hor_path_list[1:]
self.act_hor_path_list.sort(cmp=lambda a,b:cmp(a.points[0].x,b.points[0].x))
self.act_hor_path_list.sort(cmp=lambda a, b:
cmp(a.points[0].x, b.points[0].x))
scanline = []
i = 0
while i < len(self.act_hor_path_list):
s = self.act_hor_path_list[i]
if DEBUG_POLYGONEXTRACTOR2: print "s =",s
if DEBUG_POLYGONEXTRACTOR2:
print "s =", s
scanline.append(s.points[0])
if s.points[0].x <= next_x:
if len(s.points) <= 1:
if DEBUG_POLYGONEXTRACTOR2: print "remove list"
self.act_hor_path_list = self.act_hor_path_list[:i]+self.act_hor_path_list[i+1:]
if DEBUG_POLYGONEXTRACTOR2:
print "remove list"
self.act_hor_path_list = self.act_hor_path_list[:i] \
+ self.act_hor_path_list[i+1:]
# TODO: the variable "hor_list_removed" is never used.
# Any idea?
hor_list_removed = True
continue
else:
if DEBUG_POLYGONEXTRACTOR2: print "remove point", s.points[0]
if DEBUG_POLYGONEXTRACTOR2:
print "remove point", s.points[0]
s.points = s.points[1:]
if len(s.points)> 0 and s.points[0].x == next_x:
# TODO: the variable "repeat" is never used.
# Any idea?
repeat = True
i += 1
self.act_hor_path_list.sort(cmp=lambda a,b:cmp(a.points[0].x,b.points[0].x))
if len(scanline)==0:
self.act_hor_path_list.sort(cmp=lambda a, b:
cmp(a.points[0].x, b.points[0].x))
if len(scanline) == 0:
return
scanline.sort(cmp=lambda a,b:cmp(a.y,b.y))
if DEBUG_POLYGONEXTRACTOR2:
scanline.sort(cmp=lambda a, b: cmp(a.y, b.y))
if DEBUG_POLYGONEXTRACTOR2:
print "scanline' =", scanline
print "ver_hor_path_list =", self.ver_hor_path_list
print "act_hor_path_list =", self.act_hor_path_list
......
src/pycam/Geometry/Triangle.py
View file @ 107ff4dc
......@@ -21,38 +21,42 @@ You should have received a copy of the GNU General Public License
along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
"""
from Point import *
from Plane import *
from utils import *
from Line import *
from pycam.Geometry.Point import Point
from pycam.Geometry.Plane import Plane
#from pycam.Geometry.utils import *
from pycam.Geometry.Line import Line
try:
import OpenGL.GL as GL
import OpenGL.GLU as GLU
import OpenGL.GLUT as GLUT
GL_enabled = True
except:
except ImportError:
GL_enabled = False
class Triangle:
id = 0
# points are expected to be in ClockWise order
def __init__(self, p1=None, p2=None, p3=None, e1=None, e2=None, e3=None, n=None):
def __init__(self, p1=None, p2=None, p3=None, e1=None, e2=None, e3=None,
n=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)
self.e1 = Line(p1, p2)
else:
self.e1 = e1
if (not e2) and p2 and p3:
self.e2 = Line(p2,p3)
self.e2 = Line(p2, p3)
else:
self.e2 = e2
if (not e3) and p3 and p1:
self.e3 = Line(p3,p1)
self.e3 = Line(p3, p1)
else:
self.e3 = e3
self._normal = n
......@@ -69,7 +73,7 @@ class Triangle:
self._plane = None
def __repr__(self):
return "Triangle%d<%s,%s,%s>" % (self.id,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
......@@ -96,7 +100,7 @@ class Triangle:
if False and hasattr(self, "_middle"): # display bounding sphere
GL.glPushMatrix()
GL.glTranslate(self._middle.x, self._middle.y, self._middle.z)
if not hasattr(self,"_sphere"):
if not hasattr(self, "_sphere"):
self._sphere = GLU.gluNewQuadric()
GLU.gluSphere(self._sphere, self._radius, 10, 10)
GL.glPopMatrix()
......@@ -104,49 +108,51 @@ class Triangle:
GL.glPushMatrix()
cc = GL.glGetFloatv(GL.GL_CURRENT_COLOR)
c = self.center()
GL.glTranslate(c.x,c.y,c.z)
p12=self.p1.add(self.p2).mul(0.5)
p3_12=self.p3.sub(p12).normalize()
p2_1=self.p1.sub(self.p2).normalize()
pn=p2_1.cross(p3_12)
GL.glMultMatrixf((p2_1.x, p2_1.y, p2_1.z, 0, p3_12.x, p3_12.y, p3_12.z, 0, pn.x, pn.y, pn.z, 0, 0,0,0,1))
GL.glTranslate(c.x, c.y, c.z)
p12 = self.p1.add(self.p2).mul(0.5)
p3_12 = self.p3.sub(p12).normalize()
p2_1 = self.p1.sub(self.p2).normalize()
pn = p2_1.cross(p3_12)
GL.glMultMatrixf((p2_1.x, p2_1.y, p2_1.z, 0, p3_12.x, p3_12.y,
p3_12.z, 0, pn.x, pn.y, pn.z, 0, 0, 0, 0, 1))
n = self.normal().mul(0.01)
GL.glTranslatef(n.x,n.y,n.z)
GL.glScalef(0.003,0.003,0.003)
GL.glTranslatef(n.x, n.y, n.z)
GL.glScalef(0.003, 0.003, 0.003)
w = 0
for ch in str(self.id):
w += GLUT.glutStrokeWidth(GLUT.GLUT_STROKE_ROMAN, ord(ch))
GL.glTranslate(-w/2,0,0)
GL.glColor4f(1,1,1,0)
GL.glTranslate(-w/2, 0, 0)
GL.glColor4f(1, 1, 1, 0)
for ch in str(self.id):
GLUT.glutStrokeCharacter(GLUT.GLUT_STROKE_ROMAN, ord(ch))
GL.glPopMatrix()
GL.glColor4f(cc[0],cc[1],cc[2],cc[3])
GL.glColor4f(cc[0], cc[1], cc[2], cc[3])
if False: # draw point id on triangle face
cc = GL.glGetFloatv(GL.GL_CURRENT_COLOR)
c = self.center()
p12=self.p1.add(self.p2).mul(0.5)
p3_12=self.p3.sub(p12).normalize()
p2_1=self.p1.sub(self.p2).normalize()
pn=p2_1.cross(p3_12)
p12 = self.p1.add(self.p2).mul(0.5)
p3_12 = self.p3.sub(p12).normalize()
p2_1 = self.p1.sub(self.p2).normalize()
pn = p2_1.cross(p3_12)
n = self.normal().mul(0.01)
for p in (self.p1,self.p2,self.p3):
for p in (self.p1, self.p2, self.p3):
GL.glPushMatrix()
pp = p.sub(p.sub(c).mul(0.3))
GL.glTranslate(pp.x,pp.y,pp.z)
GL.glMultMatrixf((p2_1.x, p2_1.y, p2_1.z, 0, p3_12.x, p3_12.y, p3_12.z, 0, pn.x, pn.y, pn.z, 0, 0,0,0,1))
GL.glTranslatef(n.x,n.y,n.z)
GL.glScalef(0.001,0.001,0.001)
GL.glTranslate(pp.x, pp.y, pp.z)
GL.glMultMatrixf((p2_1.x, p2_1.y, p2_1.z, 0, p3_12.x, p3_12.y,
p3_12.z, 0, pn.x, pn.y, pn.z, 0, 0, 0, 0, 1))
GL.glTranslatef(n.x, n.y, n.z)
GL.glScalef(0.001, 0.001, 0.001)
w = 0
for ch in str(p.id):
w += GLUT.glutStrokeWidth(GLUT.GLUT_STROKE_ROMAN, ord(ch))
GL.glTranslate(-w/2,0,0)
GL.glColor4f(0.5,1,0.5,0)
GL.glTranslate(-w/2, 0, 0)
GL.glColor4f(0.5, 1, 0.5, 0)
for ch in str(p.id):
GLUT.glutStrokeCharacter(GLUT.GLUT_STROKE_ROMAN, ord(ch))
GL.glPopMatrix()
GL.glColor4f(cc[0],cc[1],cc[2],cc[3])
GL.glColor4f(cc[0], cc[1], cc[2], cc[3])
def normal(self):
if self._normal is None:
......@@ -158,7 +164,7 @@ class Triangle:
def plane(self):
if self._plane is None:
self._plane=Plane(self.center(), self.normal())
self._plane = Plane(self.center(), self.normal())
return self._plane
def point_inside(self, p):
......@@ -178,8 +184,9 @@ class Triangle:
# Compute barycentric coordinates
denom = dot00 * dot11 - dot01 * dot01
# originally, "u" and "v" are multiplied with "1/denom"
# we don't do this, to avoid division by zero (for triangles that are "almost" invalid)
# Originally, "u" and "v" are multiplied with "1/denom".
# We don't do this to avoid division by zero (for triangles that are
# "almost" invalid).
u = dot11 * dot02 - dot01 * dot12
v = dot00 * dot12 - dot01 * dot02
......@@ -188,32 +195,32 @@ class Triangle:
def minx(self):
if self._minx is None:
self._minx = min3(self.p1.x, self.p2.x, self.p3.x)
self._minx = min(self.p1.x, self.p2.x, self.p3.x)
return self._minx
def miny(self):
if self._miny is None:
self._miny = min3(self.p1.y, self.p2.y, self.p3.y)
self._miny = min(self.p1.y, self.p2.y, self.p3.y)
return self._miny
def minz(self):
if self._minz is None:
self._minz = min3(self.p1.z, self.p2.z, self.p3.z)
self._minz = min(self.p1.z, self.p2.z, self.p3.z)
return self._minz
def maxx(self):
if self._maxx is None:
self._maxx = max3(self.p1.x, self.p2.x, self.p3.x)
self._maxx = max(self.p1.x, self.p2.x, self.p3.x)
return self._maxx
def maxy(self):
if self._maxy is None:
self._maxy = max3(self.p1.y, self.p2.y, self.p3.y)
self._maxy = max(self.p1.y, self.p2.y, self.p3.y)
return self._maxy
def maxz(self):
if self._maxz is None:
self._maxz = max3(self.p1.z, self.p2.z, self.p3.z)
self._maxz = max(self.p1.z, self.p2.z, self.p3.z)
return self._maxz
def center(self):
......@@ -237,18 +244,25 @@ class Triangle:
return self._radiussq
def calc_circumcircle(self):
# we can't use the cached value of "normal", since we don't want the normalized value
# We can't use the cached value of "normal", since we don't want the
# normalized value.
normal = self.p2.sub(self.p1).cross(self.p3.sub(self.p2))
denom = normal.norm()
self._radius = (self.p2.sub(self.p1).norm()*self.p3.sub(self.p2).norm()*self.p3.sub(self.p1).norm())/(2*denom)
self._radius = (self.p2.sub(self.p1).norm() \
* self.p3.sub(self.p2).norm() * self.p3.sub(self.p1).norm()) \
/ (2 * denom)
self._radiussq = self._radius*self._radius
denom2 = 2*denom*denom
alpha = self.p3.sub(self.p2).normsq()*(self.p1.sub(self.p2).dot(self.p1.sub(self.p3))) / (denom2)
beta = self.p1.sub(self.p3).normsq()*(self.p2.sub(self.p1).dot(self.p2.sub(self.p3))) / (denom2)
gamma = self.p1.sub(self.p2).normsq()*(self.p3.sub(self.p1).dot(self.p3.sub(self.p2))) / (denom2)
self._middle = Point(self.p1.x*alpha+self.p2.x*beta+self.p3.x*gamma,
self.p1.y*alpha+self.p2.y*beta+self.p3.y*gamma,
self.p1.z*alpha+self.p2.z*beta+self.p3.z*gamma)
alpha = self.p3.sub(self.p2).normsq() \
* self.p1.sub(self.p2).dot(self.p1.sub(self.p3)) / denom2
beta = self.p1.sub(self.p3).normsq() \
* self.p2.sub(self.p1).dot(self.p2.sub(self.p3)) / denom2
gamma = self.p1.sub(self.p2).normsq() \
* self.p3.sub(self.p1).dot(self.p3.sub(self.p2)) / denom2
self._middle = Point(
self.p1.x * alpha + self.p2.x * beta + self.p3.x * gamma,
self.p1.y * alpha + self.p2.y * beta + self.p3.y * gamma,
self.p1.z * alpha + self.p2.z * beta + self.p3.z * gamma)
def subdivide(self, depth):
sub = []
......@@ -258,10 +272,10 @@ class Triangle:
p4 = self.p1.add(self.p2).div(2)
p5 = self.p2.add(self.p3).div(2)
p6 = self.p3.add(self.p1).div(2)
sub += Triangle(self.p1,p4,p6).subdivide(depth-1)
sub += Triangle(p6,p5,self.p3).subdivide(depth-1)
sub += Triangle(p6,p4,p5).subdivide(depth-1)
sub += Triangle(p4,self.p2,p5).subdivide(depth-1)
sub += Triangle(self.p1, p4, p6).subdivide(depth - 1)
sub += Triangle(p6, p5, self.p3).subdivide(depth - 1)
sub += Triangle(p6, p4, p5).subdivide(depth - 1)
sub += Triangle(p4, self.p2, p5).subdivide(depth - 1)
return sub
def reset_cache(self):
......
src/pycam/Geometry/TriangleKdtree.py
View file @ 107ff4dc
......@@ -20,73 +20,73 @@ You should have received a copy of the GNU General Public License
along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
"""
import math
from pycam.Geometry.kdtree import kdtree, Node
from Point import *
from Line import *
from Triangle import *
from kdtree import *
overlaptest = True
overlaptest=True
def SearchKdtree2d(kdtree, minx, maxx, miny, maxy):
if kdtree.bucket:
def SearchKdtree2d(tree, minx, maxx, miny, maxy):
if tree.bucket:
triangles = []
for n in kdtree.nodes:
global tests, hits, overlaptest
for n in tree.nodes:
if not overlaptest:
triangles.append(n.triangle)
else:
if not (n.bound[0]>maxx
or n.bound[1]<minx
or n.bound[2]>maxy
or n.bound[3]<miny):
if not (n.bound[0] > maxx
or n.bound[1] < minx
or n.bound[2] > maxy
or n.bound[3] < miny):
triangles.append(n.triangle)
return triangles
else:
if kdtree.cutdim==0:
if maxx<kdtree.minval:
if tree.cutdim == 0:
if maxx < tree.minval:
return []
elif maxx<kdtree.cutval:
return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)
elif maxx < tree.cutval:
return SearchKdtree2d(tree.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 minx>kdtree.maxval:
return SearchKdtree2d(tree.lo, minx, maxx, miny, maxy) \
+ SearchKdtree2d(tree.hi, minx, maxx, miny, maxy)
elif tree.cutdim == 1:
if minx > tree.maxval:
return []
elif minx>kdtree.cutval:
return SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
elif minx > tree.cutval:
return SearchKdtree2d(tree.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 maxy<kdtree.minval:
return SearchKdtree2d(tree.lo, minx, maxx, miny, maxy) \
+ SearchKdtree2d(tree.hi, minx, maxx, miny, maxy)
elif tree.cutdim == 2:
if maxy < tree.minval:
return []
elif maxy<kdtree.cutval:
return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)
elif maxy < tree.cutval:
return SearchKdtree2d(tree.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:
if miny>kdtree.maxval:
return SearchKdtree2d(tree.lo, minx, maxx, miny, maxy) \
+ SearchKdtree2d(tree.hi, minx, maxx, miny, maxy)
elif tree.cutdim == 3:
if miny > tree.maxval:
return []
elif miny>kdtree.cutval:
return SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
elif miny > tree.cutval:
return SearchKdtree2d(tree.hi, minx, maxx, miny, maxy)
else:
return SearchKdtree2d(kdtree.lo, minx, maxx, miny, maxy)+SearchKdtree2d(kdtree.hi, minx, maxx, miny, maxy)
return SearchKdtree2d(tree.lo, minx, maxx, miny, maxy) \
+ SearchKdtree2d(tree.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 = 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))
n.bound.append(min(t.p1.x, t.p2.x, t.p3.x))
n.bound.append(max(t.p1.x, t.p2.x, t.p3.x))
n.bound.append(min(t.p1.y, t.p2.y, t.p3.y))
n.bound.append(max(t.p1.y, t.p2.y, t.p3.y))
nodes.append(n)
kdtree.__init__(self, nodes, cutoff, cutoff_distance)
super(TriangleKdtree, self).__init__(nodes, cutoff, cutoff_distance)
def Search(self, minx, maxx, miny, maxy):
return SearchKdtree2d(self, minx, maxx, miny, maxy)
src/pycam/Geometry/__init__.py
View file @ 107ff4dc
......@@ -25,10 +25,9 @@ __all__ = ["utils", "Line", "Model", "Path", "Plane", "Point", "Triangle",
"PolygonExtractor", "TriangleKdtree", "intersection", "kdtree",
"Matrix"]
from Point import Point
from Line import Line
from Triangle import Triangle
from Path import Path
from Plane import Plane
from utils import *
from PolygonExtractor import PolygonExtractor
from pycam.Geometry.Point import Point
from pycam.Geometry.Line import Line
from pycam.Geometry.Triangle import Triangle
from pycam.Geometry.Path import Path
from pycam.Geometry.Plane import Plane
from pycam.Geometry.PolygonExtractor import PolygonExtractor
src/pycam/Geometry/intersection.py
View file @ 107ff4dc
......@@ -22,41 +22,42 @@ along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
from math import sqrt
import pycam.Geometry
from pycam.Geometry.utils import *
#import pycam.Geometry
from pycam.Utils.polynomials import poly4_roots
from pycam.Geometry.utils import INFINITE
from pycam.Geometry.Plane import Plane
from pycam.Geometry.Line import Line
from pycam.Geometry.Point import Point
def isNear(a, b):
return abs(a-b)<1e-6
return abs(a - b) < 1e-6
def isZero(a):
return isNear(a,0)
return isNear(a, 0)
def intersect_lines(xl,zl,nxl,nzl,xm,zm,nxm,nzm):
def intersect_lines(xl, zl, nxl, nzl, xm, zm, nxm, nzm):
X = None
Z = None
# print "xl=",xl,", zl=",zl,"nxl=",nxl,", nzl=",nzl,", X=", X, ", Z=",Z,", xm=",xm,",zm=",zm, ", nxm=",nxm,",nzm=",nzm
try:
if isZero(nzl) and isZero(nzm):
pass
elif isZero(nzl) or isZero(nxl):
X = xl
Z = zm + (xm-xl)*nxm/nzm
return (X,Z)
Z = zm + (xm - xl) * nxm / nzm
return (X, Z)
elif isZero(nzm) or isZero(nxm):
X = xm
Z = zl - (xm-xl)*nxl/nzl
return (X,Z)
Z = zl - (xm - xl) * nxl / nzl
return (X, Z)
else:
X = (zl-zm+(xm*nxm/nzm-xl*nxl/nzl))/(nxm/nzm-nxl/nzl)
X = (zl - zm +(xm * nxm / nzm - xl * nxl / nzl)) \
/ (nxm / nzm - nxl / nzl)
if X and xl < X and X < xm:
Z = zl + (X-xl)*nxl/nzl
return (X,Z)
except:
Z = zl + (X -xl) * nxl / nzl
return (X, Z)
except ZeroDivisionError:
pass
return (None,None)
return (None, None)
def intersect_plane_point(plane, direction, point):
denom = plane.n.dot(direction)
......@@ -71,156 +72,160 @@ def intersect_cylinder_point(center, axis, radius, radiussq, direction, point):
n = direction.cross(axis)
n.normalize()
# distance of the point to this plane
d = n.dot(point)-n.dot(center)
if d<-radius or d>radius:
return (None,None,INFINITE)
d = n.dot(point) - n.dot(center)
if abs(d) > radius:
return (None, None, INFINITE)
# ccl is on cylinder
d2 = sqrt(radiussq-d*d)
ccl = center.add(n.mul(d)).add(direction.mul(d2))
# take plane through ccl and axis
p = Plane(ccl,direction)
p = Plane(ccl, direction)
# intersect point with plane
(ccp,l)=intersect_plane_point(p, direction, point)
return (ccp,point,-l)
(ccp, l)=intersect_plane_point(p, direction, point)
return (ccp, point, -l)
def intersect_cylinder_line(center, axis, radius, radiussq, direction, edge):
d = edge.dir()
# take a plane throught the line and along the cylinder axis (1)
n = d.cross(axis)
if n.normsq()==0:
# no contact point, but should check here if cylinder *always* intersects line...
return (None,None,INFINITE)
if n.normsq() == 0:
# no contact point, but should check here if cylinder *always*
# intersects line...
return (None, None, INFINITE)
n.normalize()
# the contact line between the cylinder and this plane (1)
# is where the surface normal is perpendicular to the plane
# so line := ccl + \lambda * axis
if n.dot(direction)<0:
if n.dot(direction) < 0:
ccl = center.sub(n.mul(radius))
else:
ccl = center.add(n.mul(radius))
# now extrude the contact line along the direction, this is a plane (2)
n2 = direction.cross(axis)
if n2.normsq()==0:
# no contact point, but should check here if cylinder *always* intersects line...
return (None,None,INFINITE)
if n2.normsq() == 0:
# no contact point, but should check here if cylinder *always*
# intersects line...
return (None, None, INFINITE)
n2.normalize()
p = Plane(ccl, n2)
# intersect this plane with the line, this gives us the contact point
(cp,l) = intersect_plane_point(p, d, edge.p1)
(cp, l) = intersect_plane_point(p, d, edge.p1)
if not cp:
return (None,None,INFINITE)
# now take a plane through the contact line and perpendicular to the direction (3)
return (None, None, INFINITE)
# now take a plane through the contact line and perpendicular to the
# direction (3)
p2 = Plane(ccl, direction)
# the intersection of this plane (3) with the line throught the contact point
# the intersection of this plane (3) with the line through the contact point
# gives us the cutter contact point
(ccp,l) = intersect_plane_point(p2, direction, cp)
(ccp, l) = intersect_plane_point(p2, direction, cp)
cp = ccp.add(direction.mul(-l))
return (ccp,cp,-l)
return (ccp, cp, -l)
def intersect_circle_plane(center, radius, direction, triangle):
# let n be the normal to the plane
n = triangle.normal()
if n.dot(direction) == 0:
return (None,None,INFINITE)
return (None, None, INFINITE)
# project onto z=0
n2 = Point(n.x,n.y,0)
n2 = Point(n.x, n.y, 0)
if n2.normsq() == 0:
(cp,d) = intersect_plane_point(triangle.plane(), direction, center)
(cp, d) = intersect_plane_point(triangle.plane(), direction, center)
ccp = cp.sub(direction.mul(d))
return (ccp,cp,d)
return (ccp, cp, d)
n2.normalize()
# the cutter contact point is on the circle, where the surface normal is n
ccp = center.add(n2.mul(-radius))
# intersect the plane with a line through the contact point
(cp,d) = intersect_plane_point(triangle.plane(), direction, ccp)
return (ccp,cp,d)
(cp, d) = intersect_plane_point(triangle.plane(), direction, ccp)
return (ccp, cp, d)
def intersect_circle_point(center, axis, radius, radiussq, direction, point):
# take a plane through the base
p = Plane(center, axis)
# intersect with line gives ccp
(ccp,l) = intersect_plane_point(p, direction, point)
(ccp, l) = intersect_plane_point(p, direction, point)
# check if inside circle
if ccp and (center.sub(ccp).normsq()<=radiussq):
return (ccp,point,-l)
return (None,None,INFINITE)
if ccp and (center.sub(ccp).normsq() <= radiussq):
return (ccp, point, -l)
return (None, None, INFINITE)
def intersect_circle_line(center, axis, radius, radiussq, direction, edge):
# make a plane by sliding the line along the direction (1)
d = edge.dir()
if d.dot(axis)==0:
if direction.dot(axis)==0:
return (None,None,INFINITE)
p = Plane(center,axis)
(p1,l) = intersect_plane_point(p,direction,edge.p1)
(p2,l) = intersect_plane_point(p,direction,edge.p2)
pc = Line(p1,p2).closest_point(center)
if d.dot(axis) == 0:
if direction.dot(axis) == 0:
return (None, None, INFINITE)
p = Plane(center, axis)
(p1, l) = intersect_plane_point(p, direction, edge.p1)
(p2, l) = intersect_plane_point(p, direction, edge.p2)
pc = Line(p1, p2).closest_point(center)
d_sq = pc.sub(center).normsq()
if d_sq>radiussq:
return (None,None,INFINITE)
a = sqrt(radiussq-d_sq)
if d_sq > radiussq:
return (None, None, INFINITE)
a = sqrt(radiussq - d_sq)
d1 = p1.sub(pc).dot(d)
d2 = p2.sub(pc).dot(d)
ccp = None
cp = None
if d1>-a and d1<a:
if abs(d1) < a:
ccp = p1
cp = p1.sub(direction.mul(l))
elif d2>-a and d2<a:
elif abs(d2) < a:
ccp = p2
cp = p2.sub(direction.mul(l))
elif (d1<-a and d2>a) or (d2<-a and d1>a):
elif (d1 <- a and d2 > a) or (d2 < -a and d1 > a):
ccp = pc
cp = pc.sub(direction.mul(l))
return (ccp,cp,-l)
return (ccp, cp, -l)
n = d.cross(direction)
if n.normsq()==0:
# no contact point, but should check here if circle *always* intersects line...
return (None,None,INFINITE)
if n.normsq() == 0:
# no contact point, but should check here if circle *always* intersects
# line...
return (None, None, INFINITE)
n.normalize()
# take a plane through the base
p = Plane(center, axis)
# intersect base with line
(lp,l) = intersect_plane_point(p, d, edge.p1)
(lp, l) = intersect_plane_point(p, d, edge.p1)
if not lp:
return (None,None,INFINITE)
return (None, None, INFINITE)
# intersection of 2 planes: lp + \lambda v
v = axis.cross(n)
if v.normsq()==0:
return (None,None,INFINITE)
if v.normsq() == 0:
return (None, None, INFINITE)
v.normalize()
# take plane through intersection line and parallel to axis
n2 = v.cross(axis)
if n2.normsq()==0:
return (None,None,INFINITE)
if n2.normsq() == 0:
return (None, None, INFINITE)
n2.normalize()
# distance from center to this plane
dist = n2.dot(center)-n2.dot(lp)
distsq = dist*dist
if distsq>radiussq:
return (None,None,INFINITE)
dist = n2.dot(center) - n2.dot(lp)
distsq = dist * dist
if distsq > radiussq:
return (None, None, INFINITE)
# must be on circle
dist2 = sqrt(radiussq-distsq)
if d.dot(axis)<0:
dist2 = sqrt(radiussq - distsq)
if d.dot(axis) < 0:
dist2 = -dist2
ccp = center.sub(n2.mul(dist)).sub(v.mul(dist2))
p = Plane(edge.p1,d.cross(direction).cross(d))
(cp,l) = intersect_plane_point(p,direction,ccp)
return (ccp,cp,l)
p = Plane(edge.p1, d.cross(direction).cross(d))
(cp, l) = intersect_plane_point(p, direction, ccp)
return (ccp, cp, l)
def intersect_sphere_plane(center, radius, direction, triangle):
# let n be the normal to the plane
n = triangle.normal()
if n.dot(direction) == 0:
return (None,None,INFINITE)
return (None, None, INFINITE)
# the cutter contact point is on the sphere, where the surface normal is n
if n.dot(direction)<0:
if n.dot(direction) < 0:
ccp = center.sub(n.mul(radius))
else:
ccp = center.add(n.mul(radius))
# intersect the plane with a line through the contact point
(cp,d) = intersect_plane_point(triangle.plane(), direction, ccp)
return (ccp,cp,d)
(cp, d) = intersect_plane_point(triangle.plane(), direction, ccp)
return (ccp, cp, d)
def intersect_sphere_point(center, radius, radiussq, direction, point):
# line equation
......@@ -230,32 +235,33 @@ def intersect_sphere_point(center, radius, radiussq, direction, point):
# (1) in (2) gives a quadratic in \lambda
p0_x0 = center.sub(point)
a = direction.normsq()
b = 2*p0_x0.dot(direction)
b = 2 * p0_x0.dot(direction)
c = p0_x0.normsq() - radiussq
d = b*b-4*a*c
if d<0:
return (None,None,INFINITE)
if a<0:
l = (-b+sqrt(d))/(2*a)
d = b * b - 4 * a * c
if d < 0:
return (None, None, INFINITE)
if a < 0:
l = (-b + sqrt(d)) / (2 * a)
else:
l = (-b-sqrt(d))/(2*a)
l = (-b - sqrt(d)) / (2 * a)
# cutter contact point
ccp = point.add(direction.mul(-l))
return (ccp,point,l)
return (ccp, point, l)
def intersect_sphere_line(center, radius, radiussq, direction, edge):
# make a plane by sliding the line along the direction (1)
d = edge.dir()
n = d.cross(direction)
if n.normsq()==0:
# no contact point, but should check here if sphere *always* intersects line...
return (None,None,INFINITE)
if n.normsq() == 0:
# no contact point, but should check here if sphere *always* intersects
# line...
return (None, None, INFINITE)
n.normalize()
# calculate the distance from the sphere center to the plane
dist = -(center.dot(n)-edge.p1.dot(n))
if abs(dist)>radius:
return (None,None,INFINITE)
dist = - center.dot(n) - edge.p1.dot(n)
if abs(dist) > radius:
return (None, None, INFINITE)
# this gives us the intersection circle on the sphere
# now take a plane through the edge and perpendicular to the direction (2)
......@@ -266,96 +272,93 @@ def intersect_sphere_line(center, radius, radiussq, direction, edge):
n2.normalize()
# the contact point is on a big circle through the sphere...
dist2 = sqrt(radiussq-dist*dist)
dist2 = sqrt(radiussq - dist * dist)
# ... and it's on the plane (1)
ccp = center.add(n.mul(dist)).add(n2.mul(dist2))
# now intersect a line through this point with the plane (2)
p = Plane(edge.p1, n2)
(cp,l) = intersect_plane_point(p, direction, ccp)
return (ccp,cp,l)
(cp, l) = intersect_plane_point(p, direction, ccp)
return (ccp, cp, l)
def intersect_torus_plane(center, axis, majorradius, minorradius, direction, triangle):
def intersect_torus_plane(center, axis, majorradius, minorradius, direction,
triangle):
# take normal to the plane
n = triangle.normal()
if n.dot(direction)==0:
return (None,None,INFINITE)
if n.dot(axis)==1:
return (None,None,INFINITE)
if n.dot(direction) == 0:
return (None, None, INFINITE)
if n.dot(axis) == 1:
return (None, None, INFINITE)
# find place on torus where surface normal is n
b = n.mul(-1)
z = axis
a = b.sub(z.mul(z.dot(b)))
a_sq = a.normsq()
if a_sq<=0:
return (None,None,INFINITE)
if a_sq <= 0:
return (None, None, INFINITE)
a = a.div(sqrt(a_sq))
ccp = center.add(a.mul(majorradius)).add(b.mul(minorradius))
# find intersection with plane
p = triangle.plane()
(cp,l) = intersect_plane_point(p, direction, ccp)
return (ccp,cp,l)
(cp, l) = intersect_plane_point(p, direction, ccp)
return (ccp, cp, l)
def intersect_torus_point(center, axis, majorradius, minorradius, majorradiussq, minorradiussq, direction, point):
def intersect_torus_point(center, axis, majorradius, minorradius, majorradiussq,
minorradiussq, direction, point):
dist = 0
if direction.x==0 and direction.y==0: # drop
minlsq = sqr(majorradius-minorradius)
maxlsq = sqr(majorradius+minorradius)
l_sq = sqr(point.x-center.x)+sqr(point.y-center.y)
if (l_sq<minlsq) or (l_sq>maxlsq):
return (None,None,INFINITE)
if (direction.x == 0) and (direction.y == 0):
# drop
minlsq = (majorradius - minorradius) ** 2
maxlsq = (majorradius + minorradius) ** 2
l_sq = (point.x-center.x) ** 2 + (point.y - center.y) ** 2
if (l_sq < minlsq) or (l_sq > maxlsq):
return (None, None, INFINITE)
l = sqrt(l_sq)
z_sq = minorradiussq - sqr(majorradius - l)
z_sq = minorradiussq - (majorradius - l) ** 2
if z_sq < 0:
return (None,None,INFINITE)
return (None, None, INFINITE)
z = sqrt(z_sq)
ccp=Point(point.x,point.y,center.z-z)
dist = ccp.z-point.z
elif direction.z==0: # push
ccp = Point(point.x, point.y, center.z - z)
dist = ccp.z - point.z
elif direction.z == 0:
# push
z = point.z - center.z
if z<-minorradius or z>minorradius:
return (None,None,INFINITE)
l = majorradius + sqrt(minorradiussq - sqr(z))
if abs(z) > minorradius:
return (None, None, INFINITE)
l = majorradius + sqrt(minorradiussq - z * z)
n = axis.cross(direction)
d = n.dot(point)-n.dot(center)
if d<-l or d>l:
return (None,None,INFINITE)
a = sqrt(l*l-d*d)
d = n.dot(point) - n.dot(center)
if abs(d) > l:
return (None, None, INFINITE)
a = sqrt(l * l - d * d)
ccp = center.add(n.mul(d).add(direction.mul(a)))
ccp.z = point.z
dist = point.sub(ccp).dot(direction)
else: # general case
else:
# general case
x = point.sub(center)
v = direction.mul(-1)
x_x = x.dot(x)
x_v = x.dot(v)
x1 = Point(x.x,x.y,0)
v1 = Point(v.x,v.y,0)
x1 = Point(x.x, x.y, 0)
v1 = Point(v.x, v.y, 0)
x1_x1 = x1.dot(x1)
x1_v1 = x1.dot(v1)
v1_v1 = v1.dot(v1)
R2 = majorradiussq
r2 = minorradiussq
a = 1.0
b = 4*x_v
c = 2*(x_x)+4*sqr(x_v)+2*(R2-r2)-4*R2*v1_v1
d = 4*x_x*x_v+4*x_v*(R2-r2)-8*R2*(x1_v1)
e = sqr(x_x)+2*x_x*(R2-r2)+sqr(R2-r2)-4*R2*x1_x1
r = poly4_roots(a,b,c,d,e)
b = 4 * x_v
c = 2 * (x_x + 2 * x_v ** 2 + (R2 - r2) - 2 * R2 * v1_v1)
d = 4 * (x_x * x_v + x_v * (R2 - r2) - 2 * R2 * x1_v1)
e = (x_x) ** 2 + 2 * x_x * (R2 - r2) + (R2 - r2) ** 2 - 4 * R2 * x1_x1
r = poly4_roots(a, b, c, d, e)
if not r:
return (None,None,INFINITE)
elif len(r)==1:
l = r[0]
elif len(r)==2:
l = min(r[0],r[1])
elif len(r)==3:
l = min(min(r[0],r[1]),r[2])
elif len(r)==4:
l = min(min(r[0],r[1]),min(r[2],r[3]))
return (None, None, INFINITE)
else:
return (None,None,INFINITE)
l = min(r)
ccp = point.add(direction.mul(-l))
dist = l
return (ccp,point,dist)
return (ccp, point, dist)
src/pycam/Geometry/kdtree.py
View file @ 107ff4dc
......@@ -21,20 +21,18 @@ You should have received a copy of the GNU General Public License
along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
"""
import math
import pycam
try:
import OpenGL.GL as GL
GL_enabled = True
except:
except ImportError:
GL_enabled = False
class Node:
def __repr__(self):
s = "";
for i in range(0,len(self.bound)):
s += "%g : " % (self.bound[i])
s = ""
for bound in self.bound:
s += "%g : " % bound
return s
def find_max_spread(nodes):
......@@ -51,7 +49,7 @@ def find_max_spread(nodes):
maxval[j] = max(maxval[j], n.bound[j])
maxspreaddim = 0
maxspread = maxval[0]-minval[0]
for i in range(1,numdim):
for i in range(1, numdim):
spread = maxval[i]-minval[i]
if spread > maxspread:
maxspread = spread
......@@ -59,88 +57,13 @@ def find_max_spread(nodes):
return (maxspreaddim, maxspread)
class kdtree:
class kdtree(object):
id = 0
def __repr__(self):
if self.bucket:
if True:
return "(#%d)" % (len(self.nodes))
else:
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 not GL_enabled:
return
if self.bucket:
GL.glBegin(GL.GL_LINES)
GL.glVertex3d(minx,miny,minz)
GL.glVertex3d(minx,miny,maxz)
GL.glVertex3d(minx,maxy,minz)
GL.glVertex3d(minx,maxy,maxz)
GL.glVertex3d(maxx,miny,minz)
GL.glVertex3d(maxx,miny,maxz)
GL.glVertex3d(maxx,maxy,minz)
GL.glVertex3d(maxx,maxy,maxz)
GL.glVertex3d(minx,miny,minz)
GL.glVertex3d(maxx,miny,minz)
GL.glVertex3d(minx,maxy,minz)
GL.glVertex3d(maxx,maxy,minz)
GL.glVertex3d(minx,miny,maxz)
GL.glVertex3d(maxx,miny,maxz)
GL.glVertex3d(minx,maxy,maxz)
GL.glVertex3d(maxx,maxy,maxz)
GL.glVertex3d(minx,miny,minz)
GL.glVertex3d(minx,maxy,minz)
GL.glVertex3d(maxx,miny,minz)
GL.glVertex3d(maxx,maxy,minz)
GL.glVertex3d(minx,miny,maxz)
GL.glVertex3d(minx,maxy,maxz)
GL.glVertex3d(maxx,miny,maxz)
GL.glVertex3d(maxx,maxy,maxz)
GL.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 = kdtree.id
self.bucket = False
if nodes and len(nodes)>0:
if nodes and len(nodes) > 0:
self.dim = len(nodes[0].bound)
else:
self.dim = 0
......@@ -148,77 +71,155 @@ class kdtree:
self.cutoff = cutoff
self.cutoff_distance = cutoff_distance
if (len(nodes)<=self.cutoff):
if len(nodes) <= self.cutoff:
self.bucket = True
self.nodes = nodes
else:
(cutdim,spread) = find_max_spread(nodes)
(cutdim, spread) = find_max_spread(nodes)
if spread <= self.cutoff_distance:
self.bucket = True
self.nodes = nodes
else:
self.bucket = False
self.cutdim = cutdim
nodes.sort(cmp=lambda x,y:cmp(x.bound[cutdim],y.bound[cutdim]))
median = len(nodes)/2
nodes.sort(cmp=lambda x, y:
cmp(x.bound[cutdim], y.bound[cutdim]))
median = len(nodes) / 2
self.minval = nodes[0].bound[cutdim]
self.maxval = nodes[-1].bound[cutdim]
self.cutval = nodes[median].bound[cutdim]
self.lo = kdtree(nodes[0:median], cutoff, cutoff_distance)
self.hi = kdtree(nodes[median:], cutoff, cutoff_distance)
def __repr__(self):
if self.bucket:
if True:
return "(#%d)" % (len(self.nodes))
else:
s = "("
for n in self.nodes:
if len(s) > 1:
s += ", %s)" % str(n.p.id)
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 not GL_enabled:
return
if self.bucket:
GL.glBegin(GL.GL_LINES)
GL.glVertex3d(minx, miny, minz)
GL.glVertex3d(minx, miny, maxz)
GL.glVertex3d(minx, maxy, minz)
GL.glVertex3d(minx, maxy, maxz)
GL.glVertex3d(maxx, miny, minz)
GL.glVertex3d(maxx, miny, maxz)
GL.glVertex3d(maxx, maxy, minz)
GL.glVertex3d(maxx, maxy, maxz)
GL.glVertex3d(minx, miny, minz)
GL.glVertex3d(maxx, miny, minz)
GL.glVertex3d(minx, maxy, minz)
GL.glVertex3d(maxx, maxy, minz)
GL.glVertex3d(minx, miny, maxz)
GL.glVertex3d(maxx, miny, maxz)
GL.glVertex3d(minx, maxy, maxz)
GL.glVertex3d(maxx, maxy, maxz)
GL.glVertex3d(minx, miny, minz)
GL.glVertex3d(minx, maxy, minz)
GL.glVertex3d(maxx, miny, minz)
GL.glVertex3d(maxx, maxy, minz)
GL.glVertex3d(minx, miny, maxz)
GL.glVertex3d(minx, maxy, maxz)
GL.glVertex3d(maxx, miny, maxz)
GL.glVertex3d(maxx, maxy, maxz)
GL.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 dist(self, n1, n2):
dist = 0
for i in range(0,len(n1.bound)):
d = n1.bound[i]-n2.bound[i]
for i in range(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:
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:
d = self.dist(n, node)
if d < bestdist:
best = n
bestdist = d
return (best,bestdist)
return (best, bestdist)
else:
if node.bound[self.cutdim] <= self.cutval:
(best,bestdist) = self.lo.nearest_neighbor(node, dist)
(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)
(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)
(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)
(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:
if self.dim == 0:
self.dim = len(node.bound)
if self.bucket:
self.nodes.append(node)
if (len(self.nodes)>self.cutoff):
if len(self.nodes) > self.cutoff:
self.bucket = False
(cutdim,spread) = find_max_spread(self.nodes)
(cutdim, spread) = find_max_spread(self.nodes)
self.cutdim = cutdim
self.nodes.sort(cmp=lambda x,y:cmp(x.bound[cutdim],y.bound[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)
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)
......
src/pycam/Geometry/utils.py
View file @ 107ff4dc
......@@ -26,24 +26,3 @@ epsilon = 0.0001
def sqr(x):
return x*x
def min3(x,y,z):
if x<y:
xy = x
else:
xy = y
if xy<z:
return xy
else:
return z
def max3(x,y,z):
if x>y:
xy = x
else:
xy = y
if xy>z:
return xy
else:
return z
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment