In [1]:
__author__ = "Víctor Adolfo Gallego Alcalá, Ana María Martínez Gómez"
%matplotlib qt

In [2]:
import matplotlib.pyplot as plt
from __future__ import division
import numpy as np

In [3]:
class Point:
    
    def __init__(self, x=None, y=None):
        """ Initialize the point by giving it its x and y coordinates
        """
        self.x = x
        self.y = y
        
    def __repr__(self):
        """ Returns the representation of the point
        
        Returns
        -------
        The string which represents the point: '(x,y)'
        
        """
        return 'Point(x=%s, y=%s)' % (self.x, self.y)
    
    def get_array(self):
        """ Returns the array with represents the point
        
        Returns
        -------
        An array [x,y]
        
        """
        return np.asarray([self.x, self.y])
    
    def isLeft(self, p2, p3):
        """ This function computes the sign of the determinant of
        
         self.x self.y 1
         p2.x p2.y 1
         p3.x p3.y 1
         
        which represents the sign of the signed area of the triangle selfp2p3
        
        Parameters
        ----------
        self: first point
        p2: second point
        p3: third point
        
        Returns
        -------
        The sign of the signed area of the triangle selfp2p3
        
        """  
        return np.sign(-p2.x*self.y + p3.x*self.y + self.x*p2.y - p3.x*p2.y - self.x*p3.y + p2.x*p3.y)

def segmentIntersection(a1, a2, b1, b2):
    """ Calculate the intersection of two segments
    
    Parameters
    ----------
    a1: first end of the first segment
    a2: second end of the first segment
    b1: first end of the second segment
    b2: second end of the second segment
        
    Returns
    -------
    The intersection Point or nothing if the two segments given do not intersect

    """
    # We embedded the points in P^2
    u = np.cross([a1.x, a1.y, 1], [a2.x, a2.y, 1])
    v = np.cross([b1.x, b1.y, 1], [b2.x, b2.y, 1])
    
    # w is the intersection in P^2 of the two lines
    w = np.cross(u, v)
    
    c1 = a1.isLeft(b1, b2) * a2.isLeft(b2, b1)
    c2 = b1.isLeft(a1, a2) * b2.isLeft(a2, a1)

    # iff a1 and a2 are in different side of the second segment and 
    # b1 and b2 are in different side of the first segment, the segments intersect
    if c1 >= 0 and c2 >= 0:
        p = w[:2]/w[2]
        return Point(p[0], p[1])
    else:
        return
    
def intersect(P, Q, eps):
    """ Calculate the intersection of two Bézier curves
    
    Parameters
    ----------
    P: Array with the control points of the first Bézier curve
    Q: Array with the control points of the second Bézier curve
    eps: tolerance
        
    Returns
    -------
    An Array with the intersection Points, as there can be more than one. 
    It can also be empty is the curves do not intersect

    """
    m = P.shape[1] - 1 # number of control points in P
    n = Q.shape[1] - 1 # number of control points in Q

    # We order de control points by the x coordinate
    Pmin, Pmax = np.min(P, axis=1), np.max(P, axis=1)
    Qmin, Qmax = np.min(Q, axis=1), np.max(Q, axis=1)
    
    #We check if it is still possible to calculate the intersection
    if not (Pmin[0] > Qmax[0] or Pmin[1] >  Qmax[1] or Qmin[0] > Pmax[0] or Qmin[1] > Pmax[1]) :
        
        if m*(m-1)*np.max(np.linalg.norm(np.diff(P, n=2, axis=1), axis=0)) > eps:
            
            # We compute the compound Bézier polygon of P
            X = np.copy(P)
            T = np.atleast_2d(X)
            P_first = np.zeros([2, m + 1])
            P_first[:, 0] = P[:, 0]
            for i in xrange(m):
                T[:, 0:(m - i)] = 0.5 * (T[:, 0:(m - i)] + T[:, 1:(m - i) + 1])
                P_first[:, i + 1] = T[:, 0]
            int1 = intersect(P_first, Q, eps)
            int2 = intersect(T, Q, eps)
            return int1 + int2
        
        else:
            
            if n*(n-1)*np.max(np.linalg.norm(np.diff(Q, n=2, axis=1), axis=0)) > eps:
        
                #We compute the compound Bézier polygon of Q
                X = np.copy(Q)
                T = np.atleast_2d(X)
                Q_first = np.zeros([2, n + 1])
                Q_first[:, 0] = Q[:, 0]
                for i in xrange(n):
                    T[:, 0:(n - i)] = 0.5 * (T[:, 0:(n - i)] + T[:, 1:(n - i) + 1])
                    Q_first[:, i + 1] = T[:, 0]
                int1 = intersect(P, Q_first, eps)
                int2 = intersect(P, T, eps)
                return int1 + int2
            
            else:
                return [segmentIntersection(Point(P[0,0], P[1,0]), Point(P[0,-1], P[1,-1]), Point(Q[0,0], Q[1,0]), Point(Q[0,-1], Q[1,-1]))]
                               
                
    else: 
        return []
    
def parsePoints(L):
    """ Eliminate the Nones of an Array
    
    Parameters
    ----------
    L: Array
        
    Returns
    -------
    An Array which has the same elements that L has except the Nones

    """
    return [l for l in L if l != None]

In [4]:
# Ploting methods - The ones used in the bezier project

def bezier_subdivision(P, k, epsilon, lines=False):
    """Wrapper method of bezier_subdivision_aux. 
    
    Parameters
    ----------
    P : (n+1)xdim array of float. The control points of the curve.
    k : integer. Number of iterations.
    epsilon : float. Tolerance threshold.
    lines : boolean. Iff True, just returns the extreme points.
    
    Returns
    -------
    Nxdim array of float. The curve points. N is variable.
    """
    P = np.asarray(P, dtype=np.float64).T
    solution = bezier_subdivision_aux(P, k, epsilon, lines)
    return solution

def bezier_subdivision_aux(P, k, epsilon, lines=False):
    """It approximates the curve by performing sucesive subdivisions
    in the domain's interval. See Section 3.5 from Prautzsch H.,et al.
    Bézier and B-Spline Techniques for more details.
    
    Parameters
    ----------
    P : dimx(n+1) array of float. The control points of the curve.
    k : integer. Number of iterations.
    epsilon : float. Tolerance threshold.
    lines : boolean. Iff True, just returns the extreme points.
    
    Returns
    -------
    Nxdim array of float. The curve points. N is variable.
    """
    dim, n = P.shape
    n -= 1
    if n == 1:
        return P.T

    delta2_b = np.diff(P, n=2, axis=1)
    norm_delta2 = np.max(np.linalg.norm(delta2_b, axis=0))

    if lines and n * (n - 1) / 8 * norm_delta2 < epsilon:
        return np.array([P[0], P[-1]])

    if k == 0 or norm_delta2 < epsilon:
        return P.T
    b_first = np.zeros([dim, n + 1])
    X = np.copy(P)
    T = np.atleast_2d(X)
    b_first[:, 0] = P[:, 0]
    for i in xrange(n):
        T[:, 0:(n - i)] = 0.5 * (T[:, 0:(n - i)] + T[:, 1:(n - i) + 1])
        b_first[:, i + 1] = T[:, 0]
    z1 = bezier_subdivision(b_first.T, k - 1, epsilon, lines)
    z2 = bezier_subdivision(T[:, :].T, k - 1, epsilon, lines)
    return np.vstack([z1[:-1, :], z2])

In [5]:
class DrawLine():
    """ Based on Antonio Valdés code for plotting and moving a line
    """
    def __init__(self):
        self.figure, self.axes = plt.subplots(figsize=(10, 10))
        self.axes.set_xlim(-10, 10)
        self.axes.set_ylim(-10, 10)
        self.touched_circle = None
        self.circles = [[],[]]
        self.intersections = []
        self.line = [None, None]

        #We define the actions to do for every event
        self.cid_press = self.figure.canvas.mpl_connect('button_press_event', self.click_event)
        self.cid_move = self.figure.canvas.mpl_connect('motion_notify_event', self.motion_event)
        self.cid_release = self.figure.canvas.mpl_connect('button_release_event', self.release_event)
        
    def click_event(self, event):
        """ Action to do when clicking. It paints red circles when clicking with the right 
        button and blue ones otherwhise.
    
        Parameters
        ----------
        self: Object of the class DrawLine
        event: Event
        """
        flag = 0
        if event.button == 3: # right click
            flag = 1
        self.initial_event = event
        
        # Just do nothing if a control point is touched
        for c in self.axes.artists:
            if c.contains(event)[0]:
                self.touched_circle = c
                self.x0, self.y0 = c.center
                return
            
        c = plt.Circle((event.xdata, event.ydata), radius=0.2, color=(flag, 0, 1-flag))
        self.circles[flag].append(c)
        if not self.line[flag]:
            self.line[flag] = plt.Line2D(*zip(*map(lambda x: x.center, self.circles[flag])), color=(flag, 0, 1-flag))
            self.axes.add_line(self.line[flag])
        else:
            P = []
            for c in self.circles[flag]:
                P.append([c.center[0], c.center[1]])
            W = bezier_subdivision(P, 5, 0.1)
            W = W.T.tolist()
            self.line[flag].set_data(W)
            
        self.axes.add_artist(c)
        self.figure.canvas.draw()
        
    def motion_event(self, event):
        """ Action to do when moving. 
        It repaints the curve and paints the intersection between the two curves too.
    
        Parameters
        ----------
        self: Object of the class DrawLine
        event: Event
        
        """
        flag = 0
        if event.button == 3: # right click
            flag = 1
        
        # Do nothing if nothing is touched
        if self.touched_circle == None:
            return
        
        if event.inaxes == self.axes:
            dx = event.xdata - self.initial_event.xdata
            dy = event.ydata - self.initial_event.ydata
            self.touched_circle.center = self.x0 + dx, self.y0 + dy
            P = []
            for c in self.circles[flag]:
                P.append([c.center[0], c.center[1]])
            W = bezier_subdivision(P, 5, 0.1)
            W = W.T.tolist()
            self.line[flag].set_data(W)
            
            Q = []
            for c in self.circles[1-flag]:
                Q.append([c.center[0], c.center[1]])
            points = parsePoints(intersect(np.asarray(P).T, np.asarray(Q).T, 0.01))

            for p in points:
                self.intersections.append(plt.Circle((p.x, p.y), radius=0.1, color=(0, 1, 0)))
                self.axes.add_artist(self.intersections[-1])
            
            
            self.figure.canvas.draw()
        
    def release_event(self, event):
        """ When release the touched circule is set to None
    
        Parameters
        ----------
        self: Object of the class DrawLine
        event: Event
        
        """
        self.touched_circle = None
        
        
dl = DrawLine()

In [ ]: