fixed many float comparisons by including the "epsilon" inaccuracy compensation

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

src/pycam/Cutters/BaseCutter.py

@@ -100,13 +100,13 @@ class BaseCutter(object):
 
     def drop(self, triangle):
         # check bounding box collision
-        if self.minx > triangle.maxx:
+        if self.minx > triangle.maxx + epsilon:
             return None
-        if self.maxx < triangle.minx:
+        if self.maxx < triangle.minx - epsilon:
             return None
-        if self.miny > triangle.maxy:
+        if self.miny > triangle.maxy + epsilon:
             return None
-        if self.maxy < triangle.miny:
+        if self.maxy < triangle.miny - epsilon:
             return None
 
         # check bounding circle collision

src/pycam/Cutters/CylindricalCutter.py

@@ -175,7 +175,7 @@ class CylindricalCutter(BaseCutter):
         if cp:
             # check if the contact point is between the endpoints
             m = cp.sub(edge.p1).dot(edge.dir)
-            if (m < 0) or (m > edge.len):
+            if (m < 0) or (m > edge.len + epsilon):
                 return (None, INFINITE)
         return (cl, l)
 
@@ -190,7 +190,7 @@ class CylindricalCutter(BaseCutter):
 
     def intersect_cylinder_vertex(self, direction, point):
         (cl, ccp, cp, l) = self.intersect_cylinder_point(direction, point)
-        if ccp and ccp.z < self.center.z:
+        if ccp and ccp.z < self.center.z - epsilon:
             return (None, INFINITE)
         return (cl, l)
 
@@ -208,9 +208,9 @@ class CylindricalCutter(BaseCutter):
         if not ccp:
             return (None, INFINITE)
         m = cp.sub(edge.p1).dot(edge.dir)
-        if (m < 0) or (m > edge.len):
+        if (m < 0) or (m > edge.len + epsilon):
             return (None, INFINITE)
-        if ccp.z < self.center.z:
+        if ccp.z < self.center.z - epsilon:
             return (None, INFINITE)
         return (cl, l)
 

src/pycam/Cutters/SphericalCutter.py

@@ -187,7 +187,7 @@ class SphericalCutter(BaseCutter):
             # check if the contact point is between the endpoints
             d = edge.p2.sub(edge.p1)
             m = cp.sub(edge.p1).dot(d)
-            if (m < 0) or (m > d.normsq):
+            if (m < 0) or (m > d.normsq + epsilon):
                 return (None, INFINITE)
         return (cl, l)
 
@@ -202,7 +202,7 @@ class SphericalCutter(BaseCutter):
 
     def intersect_cylinder_vertex(self, direction, point):
         (cl, ccp, cp, l) = self.intersect_cylinder_point(direction, point)
-        if ccp and ccp.z < self.center.z:
+        if ccp and ccp.z < self.center.z - epsilon:
             return (None, INFINITE)
         return (cl, l)
 
@@ -220,9 +220,9 @@ class SphericalCutter(BaseCutter):
         if not ccp:
             return (None, INFINITE)
         m = cp.sub(edge.p1).dot(edge.dir)
-        if (m < 0) or (m > edge.len):
+        if (m < 0) or (m > edge.len + epsilon):
             return (None, INFINITE)
-        if ccp.z < self.center.z:
+        if ccp.z < self.center.z - epsilon:
             return (None, INFINITE)
         return (cl, l)
 

src/pycam/Cutters/ToroidalCutter.py

@@ -191,7 +191,7 @@ class ToroidalCutter(BaseCutter):
 
     def intersect_cylinder_vertex(self, direction, point):
         (cl, ccp, cp, l) = self.intersect_cylinder_point(direction, point)
-        if ccp and ccp.z < self.center.z:
+        if ccp and ccp.z < self.center.z - epsilon:
             return (None, INFINITE)
         return (cl, l)
 
@@ -206,11 +206,11 @@ class ToroidalCutter(BaseCutter):
 
     def intersect_cylinder_edge(self, direction, edge):
         (cl, ccp, cp, l) = self.intersect_cylinder_line(direction, edge)
-        if ccp and ccp.z < self.center.z:
+        if ccp and ccp.z < self.center.z - epsilon:
             return (None, INFINITE)
         if ccp:
             m = cp.sub(edge.p1).dot(edge.dir)
-            if (m < 0) or (m > edge.len):
+            if (m < 0) or (m > edge.len + epsilon):
                 return (None, INFINITE)
         return (cl, l)
 
@@ -256,7 +256,7 @@ class ToroidalCutter(BaseCutter):
         if cp:
             # check if the contact point is between the endpoints
             m = cp.sub(edge.p1).dot(edge.dir)
-            if (m < 0) or (m > edge.len):
+            if (m < 0) or (m > edge.len + epsilon):
                 return (None, INFINITE)
         return (cl, l)
 

src/pycam/Geometry/Line.py

@@ -164,14 +164,14 @@ class Line(TransformableContainer):
                 return x2, 1
             else:
                 return None, None
-        if infinite_lines or (0 <= factor <= 1):
+        if infinite_lines or (-epsilon <= factor <= 1 + epsilon):
             intersection = x1.add(a.mul(factor))
             # check if the intersection is between x3 and x4
             if infinite_lines:
                 return intersection, factor
-            elif (min(x3.x, x4.x) <= intersection.x <= max(x3.x, x4.x)) \
-                    and (min(x3.y, x4.y) <= intersection.y <= max(x3.y, x4.y)) \
-                    and (min(x3.z, x4.z) <= intersection.z <= max(x3.z, x4.z)):
+            elif (min(x3.x, x4.x) - epsilon <= intersection.x <= max(x3.x, x4.x) + epsilon) \
+                    and (min(x3.y, x4.y) - epsilon <= intersection.y <= max(x3.y, x4.y) + epsilon) \
+                    and (min(x3.z, x4.z) - epsilon <= intersection.z <= max(x3.z, x4.z) + epsilon):
                 return intersection, factor
             else:
                 # intersection outside of the length of line(x3, x4)
@@ -181,13 +181,9 @@ class Line(TransformableContainer):
             return None, None
 
     def get_cropped_line(self, minx, maxx, miny, maxy, minz, maxz):
-        if (minx <= self.minx <= self.maxx <= maxx) \
-                and (miny <= self.miny <= self.maxy <= maxy) \
-                and (minz <= self.minz <= self.maxz <= maxz):
+        if self.is_completely_inside(minx, maxx, miny, maxy, minz, maxz):
             return Line(line.p1, line.p2)
-        elif (maxx < self.minx) or (self.maxx < minx) \
-                or (maxy < self.miny) or (self.maxy < miny) \
-                or (maxz < self.minz) or (self.maxz < minz):
+        elif self.is_completely_outside(minx, maxx, miny, maxy, minz, maxz):
             return None
         else:
             # the line needs to be cropped
@@ -207,25 +203,23 @@ class Line(TransformableContainer):
                     for plane in planes]
             # remove all intersections outside the box and outside the line
             valid_intersections = [(cp, dist) for cp, dist in intersections
-                    if cp and (0 <= dist <= self.len) \
-                            and (minx <= cp.x <= maxx) \
-                            and (miny <= cp.y <= maxy) \
-                            and (minz <= cp.z <= maxz)]
+                    if cp and (-epsilon <= dist <= self.len + epsilon) and \
+                            cp.is_inside(minx, maxx, miny, maxy, minz, maxz)]
             # sort the intersections according to their distance to self.p1
             valid_intersections.sort(
                     cmp=lambda (cp1, l1), (cp2, l2): cmp(l1, l2))
             # Check if p1 is within the box - otherwise use the closest
             # intersection. The check for "valid_intersections" is necessary
             # to prevent an IndexError due to floating point inaccuracies.
-            if (minx <= self.p1.x <= maxx) and (miny <= self.p1.y <= maxy) \
-                    and (minz <= self.p1.z <= maxz) or not valid_intersections:
+            if self.p1.is_inside(minx, maxx, miny, maxy, minz, maxz) \
+                    or not valid_intersections:
                 new_p1 = self.p1
             else:
                 new_p1 = valid_intersections[0][0]
             # Check if p2 is within the box - otherwise use the intersection
             # most distant from p1.
-            if (minx <= self.p2.x <= maxx) and (miny <= self.p2.y <= maxy) \
-                    and (minz <= self.p2.z <= maxz) or not valid_intersections:
+            if self.p2.is_inside(minx, maxx, miny, maxy, minz, maxz) \
+                    or not valid_intersections:
                 new_p2 = self.p2
             else:
                 new_p2 = valid_intersections[-1][0]

src/pycam/Geometry/Matrix.py

@@ -24,7 +24,7 @@ along with PyCAM.  If not, see <http://www.gnu.org/licenses/>.
 
 
 from pycam.Geometry.Point import Point
-from pycam.Geometry.utils import sqrt, number
+from pycam.Geometry.utils import sqrt, number, epsilon
 import math
 
 
@@ -112,9 +112,9 @@ def get_rotation_matrix_from_to(v_orig, v_dest):
     cross_product = get_length(get_cross_product(v_orig, v_dest))
     arcsin = cross_product / (v_orig_length * v_dest_length)
     # prevent float inaccuracies to crash the calculation (within limits)
-    if arcsin > 1.0001:
+    if 1 < arcsin < 1 + epsilon:
         arcsin = 1.0
-    elif arcsin < 1.001:
+    elif -1 - epsilon < arcsin < -1:
         arcsin = -1.0
     rot_angle = math.asin(arcsin)
     # calculate the rotation axis

src/pycam/Geometry/Point.py

@@ -21,11 +21,11 @@ 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.Geometry.utils import sqrt, number
+from pycam.Geometry.utils import epsilon, sqrt, number
 from decimal import Decimal
 
 def _is_near(x, y):
-    return abs(x - y) < 1e-6
+    return abs(x - y) < epsilon
 
 
 class Point:
@@ -105,3 +105,12 @@ class Point:
         else:
             return Point(self.x / n, self.y / n, self.z / n)
 
+    def is_inside(self, minx=None, maxx=None, miny=None, maxy=None, minz=None,
+            maxz=None):
+        return ((minx is None) or (minx - epsilon <= self.x)) \
+                and ((maxx is None) or (self.x <= maxx + epsilon)) \
+                and ((miny is None) or (miny - epsilon <= self.y)) \
+                and ((maxy is None) or (self.y <= maxy + epsilon)) \
+                and ((minz is None) or (minz - epsilon <= self.z)) \
+                and ((maxz is None) or (self.z <= maxz + epsilon))
+

src/pycam/Geometry/Polygon.py

@@ -24,7 +24,7 @@ from pycam.Geometry.Line import Line
 from pycam.Geometry.Point import Point
 from pycam.Geometry.Plane import Plane
 from pycam.Geometry import TransformableContainer
-from pycam.Geometry.utils import number
+from pycam.Geometry.utils import number, epsilon
 import pycam.Geometry.Matrix as Matrix
 import math
 
@@ -154,9 +154,8 @@ class Polygon(TransformableContainer):
         The result is unpredictable if the point is exactly on one of the lines.
         """
         # First: check if the point is within the boundary of the polygon.
-        if not ((self.minx <= p.x <= self.maxx) \
-                and (self.miny <= p.y <= self.maxy) \
-                and (self.minz <= p.z <= self.maxz)):
+        if not p.is_inside(self.minx, self.maxx, self.miny, self.maxy,
+                self.minz, self.maxz):
             # the point is outside the rectangle boundary
             return False
         # see http://www.alienryderflex.com/polygon/
@@ -175,7 +174,7 @@ class Polygon(TransformableContainer):
                     or ((p2.y < p.y) and (p.y <= p1.y)):
                 part_y = (p.y - p1.y) / (p2.y - p1.y)
                 intersection_x = p1.x + part_y * (p2.x - p1.x)
-                if intersection_x < p.x:
+                if intersection_x < p.x + epsilon:
                     # only count intersections to the left
                     intersection_count += 1
         # odd intersection count -> inside
@@ -598,9 +597,7 @@ class Polygon(TransformableContainer):
         new_groups = []
         for line in self.get_lines():
             new_line = None
-            if (minx <= line.minx <= line.maxx <= maxx) \
-                    and (miny <= line.miny <= line.maxy <= maxy) \
-                    and (minz <= line.minz <= line.maxz <= maxz):
+            if line.is_completely_inside(minx, maxx, miny, maxy, minz, maxz):
                 new_line = line
             else:
                 cropped_line = line.get_cropped_line(minx, maxx, miny, maxy,

src/pycam/Geometry/__init__.py

@@ -25,6 +25,8 @@ __all__ = ["utils", "Line", "Model", "Path", "Plane", "Point", "Triangle",
            "PolygonExtractor", "TriangleKdtree", "intersection", "kdtree",
            "Matrix", "Polygon"]
 
+from pycam.Geometry.utils import epsilon
+
 
 class TransformableContainer(object):
     """ a base class for geometrical objects containing other elements
@@ -95,3 +97,21 @@ class TransformableContainer(object):
                 + "'TransformableContainer' but it fails to implement the " \
                 + "'reset_cache' method") % str(type(self)))
 
+    def is_completely_inside(self, minx=None, maxx=None, miny=None, maxy=None,
+            minz=None, maxz=None):
+        return ((minx is None) or (minx - epsilon <= self.minx)) \
+                and ((maxx is None) or (self.maxx <= maxx + epsilon)) \
+                and ((miny is None) or (miny - epsilon <= self.miny)) \
+                and ((maxy is None) or (self.maxy <= maxy + epsilon)) \
+                and ((minz is None) or (minz - epsilon <= self.minz)) \
+                and ((maxz is None) or (self.maxz <= maxz + epsilon))
+
+    def is_completely_outside(self, minx=None, maxx=None, miny=None, maxy=None,
+            minz=None, maxz=None):
+        return ((maxx is None) or (maxx + epsilon < self.minx)) \
+                or ((minx is None) or (self.maxx < minx - epsilon)) \
+                or ((maxy is None) or (maxy + epsilon < self.miny)) \
+                or ((miny is None) or (self.maxy < miny - epsilon)) \
+                or ((maxz is None) or (maxz + epsilon < self.minz)) \
+                or ((minz is None) or (self.maxz < minz - epsilon))
+

src/pycam/Geometry/intersection.py

@@ -23,13 +23,13 @@ along with PyCAM.  If not, see <http://www.gnu.org/licenses/>.
 
 #import pycam.Geometry
 from pycam.Utils.polynomials import poly4_roots
-from pycam.Geometry.utils import INFINITE, sqrt
+from pycam.Geometry.utils import INFINITE, sqrt, epsilon
 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) < epsilon
 
 def isZero(a):
     return isNear(a, 0)
@@ -63,7 +63,7 @@ def intersect_cylinder_point(center, axis, radius, radiussq, direction, point):
     n = direction.cross(axis).normalized()
     # distance of the point to this plane
     d = n.dot(point) - n.dot(center)
-    if abs(d) > radius:
+    if abs(d) > radius + epsilon:
         return (None, None, INFINITE)
     # ccl is on cylinder
     d2 = sqrt(radiussq-d*d)
@@ -157,13 +157,14 @@ def intersect_circle_line(center, axis, radius, radiussq, direction, edge):
         d2 = p2.sub(pc).dot(d)
         ccp = None
         cp = None
-        if abs(d1) < a:
+        if abs(d1) < a + epsilon:
             ccp = p1
             cp = p1.sub(direction.mul(l))
-        elif abs(d2) < a:
+        elif abs(d2) < a + epsilon:
             ccp = p2
             cp = p2.sub(direction.mul(l))
-        elif (d1 <- a and d2 > a) or (d2 < -a and d1 > a):
+        elif ((d1 < -a + epsilon) and (d2 > a - epsilon)) \
+                or ((d2 < -a + epsilon) and (d1 > a - epsilon)):
             ccp = pc
             cp = pc.sub(direction.mul(l))
         return (ccp, cp, -l)
@@ -300,7 +301,7 @@ def intersect_torus_point(center, axis, majorradius, minorradius, majorradiussq,
         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):
+        if (l_sq < minlsq - epsilon) or (l_sq > maxlsq + epsilon):
             return (None, None, INFINITE)
         l = sqrt(l_sq)
         z_sq = minorradiussq - (majorradius - l) ** 2
@@ -312,12 +313,12 @@ def intersect_torus_point(center, axis, majorradius, minorradius, majorradiussq,
     elif direction.z == 0:
         # push
         z = point.z - center.z
-        if abs(z) > minorradius:
+        if abs(z) > minorradius + epsilon:
             return (None, None, INFINITE)
         l = majorradius + sqrt(minorradiussq - z * z)
         n = axis.cross(direction)
         d = n.dot(point) - n.dot(center)
-        if abs(d) > l:
+        if abs(d) > l + epsilon:
             return (None, None, INFINITE)
         a = sqrt(l * l - d * d)
         ccp = center.add(n.mul(d).add(direction.mul(a)))

src/pycam/PathGenerators/__init__.py

@@ -79,7 +79,7 @@ def get_free_paths_triangles(model, cutter, p1, p2):
     for h in hits:
         if h.dir == forward:
             if count == 0:
-                if h.d >= 0:
+                if h.d >= 0 - epsilon:
                     if len(points) == 0:
                         points.append(p1)
                     points.append(h.cl)
@@ -87,7 +87,7 @@ def get_free_paths_triangles(model, cutter, p1, p2):
         else:
             count -= 1
             if count == 0:
-                if h.d <= xyz_dist:
+                if h.d <= xyz_dist + epsilon:
                     points.append(h.cl)
 
     if len(points)%2 == 1:
@@ -233,18 +233,18 @@ def get_max_height_triangles(model, cutter, x, y, minz, maxz, order=None,
             box_y_max, box_z_max)
     for t in triangles:
         cut = cutter.drop(t)
-        if cut and (cut.z > height_max or height_max is None):
+        if cut and ((height_max is None) or (cut.z > height_max)):
             height_max = cut.z
             cut_max = cut
             triangle_max = t
     # don't do a complete boundary check for the height
     # this avoids zero-cuts for models that exceed the bounding box height
-    if not cut_max or cut_max.z < minz:
+    if not cut_max or cut_max.z < minz + epsilon:
         cut_max = Point(x, y, minz)
     if last_pos["cut"] and \
             ((triangle_max and not last_pos["triangle"]) \
             or (last_pos["triangle"] and not triangle_max)):
-        if minz <= last_pos["cut"].z <= maxz:
+        if minz - epsilon <= last_pos["cut"].z <= maxz + epsilon:
             result.append(Point(last_pos["cut"].x, last_pos["cut"].y,
                     cut_max.z))
         else: