little python hacks: intervals and line segments

[Kragen Sitaker]

An interval-math ADT and a line-segment ADT.  Insomnia.

# Kragen Sitaker, 2000-12-30.
# Written pursuant to SICP data abstraction exercise.  Not intended
# for real-life use, and here's why:
# - __mul__ (and therefore __div__) is likely unnecessarily slow
# - mixed-mode arithmetic isn't really defined; this is needed at least for
#   things like 1/x, where x is an Interval
# - it's quite likely buggy; I've only used it in cute little tests
# - there are no automated tests

class Interval:
    # general python glue methods
    def __init__(self, min, max):
        self.min = min
        self.max = max
    def __repr__(self):
        return "Interval(%s, %s)" % (self.min, self.max)
    def __str__(self):
        if self.min == self.max:
            return str(self.min)
        else:
            return "(%s..%s)" % (self.min, self.max)

    # arithmetic methods
    def __add__(self, other):
        return Interval(self.min + other.min,
                        self.max + other.max)
    def __neg__(self): return Interval(-self.max, -self.min)
    def __pos__(self): return self  # silly!
    def __abs__(self):
        if self.min <= 0 and self.max >= 0:
            return Interval(0, max([abs(self.min), abs(self.max)]))
        elif self.min <= 0:
            return Interval(abs(self.max), abs(self.min))
        else:
            return Interval(abs(self.min), abs(self.max))
    def __sub__(self, other): return self + -other
    def __mul__(self, other):
        corners = [self.min * other.min,
                   self.max * other.max,
                   self.min * other.max,
                   self.max * other.min]
        return Interval(min(corners), max(corners))
    def reciprocal(self):
        if self.min <= 0 and self.max >= 0:
            raise ZeroDivisionError, "interval division"
        else:
            return Interval(1./self.max, 1./self.min)
    def __div__(self, other): return self * other.reciprocal()

def const(value): return Interval(value, value)
def make_center_width(center, width):
    return Interval(center - width, center + width)
def tolerance(center, error):
    return make_center_width(center, error * center)


And here's the line-segment ADT:

# Kragen Sitaker, 2000-12-29.
# Representation of line segments, pursuant to a SICP exercise in
# the "Abstraction Barriers" section, although I did get a little
# carried away.
# Not intended for real use:
# * there are probably bugs
# * I'm not using it except for fun
# * it's too tied to floating-point, which has disadvantages for
#   eometrical stuff

from math import sqrt

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def __repr__(self):
        return "Point(%s, %s)" % (self.x, self.y)
    def __cmp__(self, other):
        if self.x == other.x and self.y == other.y:
            return 0
        else: # try to come up with something stable
            return cmp(id(self), id(other))

class Segment:
    def __init__(self, start, end):
        self.start = start
        self.end = end
    def midpoint(self):
        return Point((self.start.x + self.end.x)/2.,
                     (self.start.y + self.end.y)/2.)
    def length(self):
        return sqrt((self.start.x - self.end.x)**2 +
                    (self.start.y - self.end.y)**2)
    def __repr__(self):
        return "[%s..%s]" % (self.start, self.end)

    # now for some fun stuff

    def vertical(self):
        return self.start.x == self.end.x
    # don't call this if the line is vertical
    def slope(self):
        return (self.end.y - self.start.y)/(self.end.x - self.start.x)
    # or this
    def yintercept(self):
        return self.start.y - self.start.x * self.slope()
    def parallel(self, other):
        if self.vertical():
            return other.vertical()
        elif other.vertical():
            return 0
        else:
            return self.slope() == other.slope()

    # return the point at which two segments would intersect if they extended
    # far enough
    def intersection(self, other):
        if self.parallel(other):
            # the segments may intersect, but we don't care
            return None
        elif self.vertical():
            return other.intersection(self)
        elif other.vertical():
            return Point(other.start.x,
                         self.yintercept() + other.start.x * self.slope())
        else:
            # m1x + b1 = m2x + b2; so
            # (m1 - m2)x + b1 - b2 == 0
            # (m1 - m2)x = b2 - b1
            # x = (b2 - b1)/(m1 - m2)
            x = ((other.yintercept() - self.yintercept()) /
                 (self.slope() - other.slope()))
            return Point(x, self.yintercept() + x * self.slope())

    def in_bbox(self, point):
        return ((point.x <= self.start.x and point.x >= self.end.x or
                 point.x <= self.end.x and point.x >= self.start.x) and
                (point.y <= self.start.y and point.y >= self.end.y or
                 point.y <= self.end.y and point.y >= self.start.y))
    # is a point collinear with this line segment?
    def on_line(self, point):
        if self.vertical():
            return point.x == self.start.x
        else:
            return (point.x * self.slope() + self.yintercept() == point.y)

    def intersects(self, other):
        if self.parallel(other):
            # they can "intersect" if they are collinear and overlap
            if not (self.in_bbox(other.start) or self.in_bbox(other.end)):
                return 0
            elif self.vertical():
                return self.start.x == other.start.x
            else:
                return self.yintercept() == other.yintercept()
        else:
            i = self.intersection(other)
            assert i is not None
            return self.in_bbox(i) and other.in_bbox(i)
        
def test():
    x = Segment(Point(1., 2.), Point(4., 6.))
    y = Segment(Point(2., 3.), Point(5., 7.))
    z = Segment(Point(1., 1.), Point(1., 3.))
    assert x.start == Point(1., 2.)
    assert x.start != x.end
    assert x.parallel(y)
    assert y.parallel(x)
    assert not x.parallel(z)
    assert not y.parallel(z)
    assert not z.parallel(x)
    assert not z.parallel(y)
    assert x.intersection(y) is None
    assert y.intersection(x) is None
    assert z.intersection(x) is not None
    assert x.intersection(z) is not None
    i = z.intersection(y)
    assert i is not None # even though they don't intersect
    assert z.on_line(i)
    assert y.on_line(i)
    assert z.in_bbox(i)
    assert not y.in_bbox(i)
    assert z.intersection(x) == x.intersection(z)
    assert x.intersection(z) == Point(1., 2.)
    xcross = Segment(Point(1., 6.), Point(4., 2.))
    assert x.length() == xcross.length()
    assert x.length() == 5.0
    assert xcross.intersection(x) is not None
    assert xcross.intersection(x) == x.intersection(xcross)
    assert x.intersection(xcross) == x.midpoint()
    assert z.intersects(x)
    assert not z.intersects(y)
    assert not z.intersects(xcross)
    assert not x.intersects(y)
    v1 = Segment(Point(1., 10.), Point(1., 20.))
    v2 = Segment(Point(1., 3.), Point(1., 20.))
    v3 = Segment(Point(2., 1.), Point(2., 3.))
    assert z.intersection(v1) is None
    assert not z.intersects(v1)
    assert z.intersects(v2)
    assert not z.intersects(v3)
    assert v1.intersects(v2)
    assert not v1.intersects(v3)

test()