Crop a collection of lines using another collection of lines

Question

First of all, here is a quick graphical description of what I intend to do. I will use Python, but feel free to use pseudo code in your answers.

I have 2 collections of 2D segments, stored as the following: [ [start_point, end_point], [...] ].

For the first step, I have to detect each segment of the blue collection which is colliding with the black collection. For this I use LeMothe's line/line intersection algorithm.

Then comes my first problem: Having a segment AB which intersects with my line in C, I don't know how to determine using code if I have to trim my segment as AC or as CB , ie: I don't know which part of my segment I need to keep and which one I need to remove.

Then for the second step, I have really no ideas how to achieve this.

Any help would be greatly appreciated, so thank you in advance!

Answer 1

The second step is trivial once you figure what to keep and what not, you just need to keep track of the segments you clipped and see where they were originally joined (e.g. assume that the segments are in order and form a connected line). On the other hand, given that your black line is in fact a line and not a polygon, in your first step, choosing what is "outside" and what is "inside" seems completely arbitrary; is it possible to close that into a polygon? Otherwise, you may need to artificially create two polygons (one for each side of the line) and then do clipping inside those polygons. You could use something like the Cyrus and Beck line clipping algorithm (see this tutorial for an overview: https://www.tutorialspoint.com/computer_graphics/viewing_and_clipping.htm)

Feel free to use any of the code below as a starting point (you have an intersect function and some useful classes). Implements Sutherland and Hodgman.

class Point2(object):
    """Structure for a 2D point"""
    def __init__(self, x=0, y=0):
        self.x = x
        self.y = y
    def __copy__(self):
        return self.__class__(self.x, self.y)
    copy = __copy__
    def __repr__(self):
        return 'Point2(%d, %d)' % (self.x, self.y)
    def __getitem__(self, key):
        return (self.x, self.y)[key]
    def __setitem__(self, key, value):
        l = [self.x, self.y]
        l[key] = value
        self.x, self.y = l
    def __eq__(self, other):
        if isinstance(other, Point2):
            return self.x == other.x and \
                   self.y == other.y
        else:
            assert hasattr(other, '__len__') and len(other) == 2
            return self.x == other[0] and \
                   self.y == other[1]
    def __ne__(self, other):
        return not self.__eq__(other)
    def __nonzero__(self):
        return self.x != 0 or self.y != 0
    def __len__(self):
        return 2


class Line2(object):
    """Structure for a 2D line"""
    def __init__(self,pt1,pt2):
        self.pt1,self.pt2=pt1,pt2
    def __repr__(self):
        return 'Line2(%s, %s)' % (self.pt1, self.pt2)

class Polygon2(object):
    def __init__(self,points):
        self.points = points
    def __repr__(self):
        return '[\n %s\n]' % '\n '.join([str(i) for i in self.points])
    def lines(self):
        lines = []
        e = self.points[-1].copy()
        for p in self.points:
            lines.append(Line2(e,p))
            e = p.copy()
        return lines
        #return [Line2(a,b) for a,b in zip(self.points,self.points[1:]+[self.points[0]])]
    def __copy__(self):
        return self.__class__(list(self.points))
    copy = __copy__


class Renderer(object):
    """Rendering algorithm implementations"""
    def __init__(self,world,img,color=1):
        self.world,self.img,self.color=world,img,color
    def transform(self,s,r,m,n):
        """Homogeneous transformation operations"""
        for i in self.world.points():
            j = Matrix3.new_translate(m, n)*Matrix3.new_rotate(r)*Matrix3.new_scale(s)*i
            i.x,i.y = j.x,j.y
    def clip(self,a,b,c,d):
        """Clipping for the world window defined by a,b,c,d"""
        self.clip_lines(a, b, c, d)
        self.clip_polygons(a, b, c, d)
    def shift(self,a,b,c,d):
        """Shift the world window"""
        for i in self.world.points():
            i.x -= a
            i.y -= b
    def clip_lines(self,a,b,c,d):
        """Clipping for lines (i.e. open polygons)"""
        clipped = []
        for i in self.world.lines:
            clipped += [self.clip_lines_cohen_sutherland(i.pt1, i.pt2, a, b, c, d)]
        self.world.lines = [i for i in clipped if i]
    def clip_polygons(self,a,b,c,d):
        """Clipping for polygons"""
        polygons = []
        for polygon in self.world.polygons:
            new_polygon = self.clip_polygon_sutherland_hodgman(polygon, a, b, c, d)
            polygons.append(new_polygon)
        self.world.polygons = polygons
    def clip_polygon_sutherland_hodgman(self,polygon,xmin,ymin,xmax,ymax):
        edges = [Line2(Point2(xmax,ymax),Point2(xmin,ymax)), #top
                 Line2(Point2(xmin,ymax),Point2(xmin,ymin)), #left
                 Line2(Point2(xmin,ymin),Point2(xmax,ymin)), #bottom
                 Line2(Point2(xmax,ymin),Point2(xmax,ymax)), #right
                 ]
        def is_inside(pt,line):
            # uses the determinant of the vectors (AB,AQ), Q(X,Y) is the query
            # left is inside
            det = (line.pt2.x-line.pt1.x)*(pt.y-line.pt1.y) - (line.pt2.y-line.pt1.y)*(pt.x-line.pt1.x)
            return det>=0
        def intersect(pt0,pt1,line):
            x1,x2,x3,x4 = pt0.x,pt1.x,line.pt1.x,line.pt2.x
            y1,y2,y3,y4 = pt0.y,pt1.y,line.pt1.y,line.pt2.y
            x = ((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4))
            y = ((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4)) / ((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4))
            return Point2(int(x),int(y))
        polygon_new = polygon.copy()
        for edge in edges:
            polygon_copy = polygon_new.copy()
            polygon_new = Polygon2([])
            s = polygon_copy.points[-1]
            for p in polygon_copy.points:
                if is_inside(s,edge) and is_inside(p,edge):
                    polygon_new.points.append(p)
                elif is_inside(s,edge) and not is_inside(p,edge):
                    polygon_new.points.append(intersect(s,p,edge))
                elif not is_inside(s,edge) and not is_inside(p,edge):
                    pass
                else:
                    polygon_new.points.append(intersect(s,p,edge))
                    polygon_new.points.append(p)
                s = p
        return polygon_new
    def clip_lines_cohen_sutherland(self,pt0,pt1,xmin,ymin,xmax,ymax):
        """Cohen-Sutherland clipping algorithm for line pt0 to pt1 and clip rectangle with diagonal from (xmin,ymin) to (xmax,ymax)."""
        TOP = 1
        BOTTOM = 2
        RIGHT = 4
        LEFT = 8
        def ComputeOutCode(pt):
            code = 0
            if pt.y > ymax: code += TOP
            elif pt.y < ymin: code += BOTTOM
            if pt.x > xmax: code += RIGHT
            elif pt.x < xmin: code += LEFT
            return code
        accept = False
        outcode0, outcode1 = ComputeOutCode(pt0), ComputeOutCode(pt1)
        while True:
            if outcode0==outcode1==0:
                accept=True
                break
            elif outcode0&outcode1:
                accept=False
                break
            else:
                #Failed both tests, so calculate the line segment to clip from an outside point to an intersection with clip edge.
                outcodeOut = outcode0 if not outcode0 == 0 else outcode1
                if TOP & outcodeOut:
                    x = pt0.x + (pt1.x - pt0.x) * (ymax - pt0.y) / (pt1.y - pt0.y)
                    y = ymax
                elif BOTTOM & outcodeOut:
                    x = pt0.x + (pt1.x - pt0.x) * (ymin - pt0.y) / (pt1.y - pt0.y)
                    y = ymin
                elif RIGHT & outcodeOut:
                    y = pt0.y + (pt1.y - pt0.y) * (xmax - pt0.x) / (pt1.x - pt0.x);
                    x = xmax;
                elif LEFT & outcodeOut:
                    y = pt0.y + (pt1.y - pt0.y) * (xmin - pt0.x) / (pt1.x - pt0.x);
                    x = xmin;
                if outcodeOut == outcode0:
                    pt0 = Point2(x,y)
                    outcode0 = ComputeOutCode(pt0)
                else:
                    pt1 = Point2(x,y)
                    outcode1 = ComputeOutCode(pt1);
        if accept:
            return Line2(pt0,pt1)
        else:
            return False