python-matplotlib 1.1.1~rc1+git20120423-0ubuntu1 / usr / share / pyshared / matplotlib / bezier.py

"""
A module providing some utility functions regarding bezier path manipulation.
"""


import numpy as np
from math import sqrt

from matplotlib.path import Path

from operator import xor
import warnings


class NonIntersectingPathException(ValueError):
    pass

# some functions

def get_intersection(cx1, cy1, cos_t1, sin_t1,
                     cx2, cy2, cos_t2, sin_t2):
    """ return a intersecting point between a line through (cx1, cy1)
    and having angle t1 and a line through (cx2, cy2) and angle t2.
    """

    # line1 => sin_t1 * (x - cx1) - cos_t1 * (y - cy1) = 0.
    # line1 => sin_t1 * x + cos_t1 * y = sin_t1*cx1 - cos_t1*cy1

    line1_rhs = sin_t1 * cx1 - cos_t1 * cy1
    line2_rhs = sin_t2 * cx2 - cos_t2 * cy2

    # rhs matrix
    a, b = sin_t1, -cos_t1
    c, d = sin_t2, -cos_t2

    ad_bc = a*d-b*c
    if ad_bc == 0.:
        raise ValueError("Given lines do not intersect")

    #rhs_inverse
    a_, b_ = d, -b
    c_, d_ = -c, a
    a_, b_, c_, d_ = [k / ad_bc for k in [a_, b_, c_, d_]]

    x = a_* line1_rhs + b_ * line2_rhs
    y = c_* line1_rhs + d_ * line2_rhs

    return x, y



def get_normal_points(cx, cy, cos_t, sin_t, length):
    """
    For a line passing through (*cx*, *cy*) and having a angle *t*,
    return locations of the two points located along its perpendicular line at the distance of *length*.
    """

    if length == 0.:
        return cx, cy, cx, cy

    cos_t1, sin_t1 = sin_t, -cos_t
    cos_t2, sin_t2 = -sin_t, cos_t

    x1, y1 = length*cos_t1 + cx, length*sin_t1 + cy
    x2, y2 = length*cos_t2 + cx, length*sin_t2 + cy

    return x1, y1, x2, y2




## BEZIER routines





# subdividing bezier curve
# http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/Bezier/bezier-sub.html

def _de_casteljau1(beta, t):
    next_beta = beta[:-1] * (1-t) + beta[1:] * t
    return next_beta

def split_de_casteljau(beta, t):
    """split a bezier segment defined by its controlpoints *beta*
    into two separate segment divided at *t* and return their control points.

    """
    beta = np.asarray(beta)
    beta_list = [beta]
    while True:
        beta = _de_casteljau1(beta, t)
        beta_list.append(beta)
        if len(beta) == 1:
            break
    left_beta = [beta[0] for beta in beta_list]
    right_beta = [beta[-1] for beta in reversed(beta_list)]

    return left_beta, right_beta







def find_bezier_t_intersecting_with_closedpath(bezier_point_at_t, inside_closedpath,
                                               t0=0., t1=1., tolerence=0.01):
    """ Find a parameter t0 and t1 of the given bezier path which
    bounds the intersecting points with a provided closed
    path(*inside_closedpath*). Search starts from *t0* and *t1* and it
    uses a simple bisecting algorithm therefore one of the end point
    must be inside the path while the orther doesn't. The search stop
    when |t0-t1| gets smaller than the given tolerence.
    value for

    - bezier_point_at_t : a function which returns x, y coordinates at *t*

    - inside_closedpath : return True if the point is insed the path

    """
    # inside_closedpath : function

    start = bezier_point_at_t(t0)
    end = bezier_point_at_t(t1)

    start_inside = inside_closedpath(start)
    end_inside = inside_closedpath(end)

    if not xor(start_inside, end_inside):
        raise NonIntersectingPathException("the segment does not seemed to intersect with the path")

    while 1:

        # return if the distance is smaller than the tolerence
        if (start[0]-end[0])**2 + (start[1]-end[1])**2 < tolerence**2:
            return t0, t1

        # calculate the middle point
        middle_t = 0.5*(t0+t1)
        middle = bezier_point_at_t(middle_t)
        middle_inside = inside_closedpath(middle)

        if xor(start_inside, middle_inside):
            t1 = middle_t
            end = middle
            end_inside = middle_inside
        else:
            t0 = middle_t
            start = middle
            start_inside = middle_inside





class BezierSegment:
    """
    A simple class of a 2-dimensional bezier segment
    """

    # Highrt order bezier lines can be supported by simplying adding
    # correcponding values.
    _binom_coeff = {1:np.array([1., 1.]),
                    2:np.array([1., 2., 1.]),
                    3:np.array([1., 3., 3., 1.])}

    def __init__(self, control_points):
        """
        *control_points* : location of contol points. It needs have a
         shpae of n * 2, where n is the order of the bezier line. 1<=
         n <= 3 is supported.
        """
        _o = len(control_points)
        self._orders = np.arange(_o)
        _coeff = BezierSegment._binom_coeff[_o - 1]

        _control_points = np.asarray(control_points)
        xx = _control_points[:,0]
        yy = _control_points[:,1]

        self._px = xx * _coeff
        self._py = yy * _coeff

    def point_at_t(self, t):
        "evaluate a point at t"
        one_minus_t_powers = np.power(1.-t, self._orders)[::-1]
        t_powers = np.power(t, self._orders)

        tt = one_minus_t_powers * t_powers
        _x = sum(tt * self._px)
        _y = sum(tt * self._py)

        return _x, _y


def split_bezier_intersecting_with_closedpath(bezier,
                                              inside_closedpath,
                                              tolerence=0.01):

    """
    bezier : control points of the bezier segment
    inside_closedpath : a function which returns true if the point is inside the path
    """

    bz = BezierSegment(bezier)
    bezier_point_at_t = bz.point_at_t

    t0, t1 = find_bezier_t_intersecting_with_closedpath(bezier_point_at_t,
                                                        inside_closedpath,
                                                        tolerence=tolerence)

    _left, _right = split_de_casteljau(bezier, (t0+t1)/2.)
    return _left, _right



def find_r_to_boundary_of_closedpath(inside_closedpath, xy,
                                     cos_t, sin_t,
                                     rmin=0., rmax=1., tolerence=0.01):
    """
    Find a radius r (centered at *xy*) between *rmin* and *rmax* at
    which it intersect with the path.

    inside_closedpath : function
    cx, cy : center
    cos_t, sin_t : cosine and sine for the angle
    rmin, rmax :
    """

    cx, cy = xy
    def _f(r):
        return cos_t*r + cx, sin_t*r + cy

    find_bezier_t_intersecting_with_closedpath(_f, inside_closedpath,
                                               t0=rmin, t1=rmax, tolerence=tolerence)



## matplotlib specific

def split_path_inout(path, inside, tolerence=0.01, reorder_inout=False):
    """ divide a path into two segment at the point where inside(x, y)
    becomes False.
    """

    path_iter = path.iter_segments()

    ctl_points, command = path_iter.next()
    begin_inside = inside(ctl_points[-2:]) # true if begin point is inside

    bezier_path = None
    ctl_points_old = ctl_points

    concat = np.concatenate

    iold=0
    i = 1

    for ctl_points, command in path_iter:
        iold=i
        i += len(ctl_points)/2
        if inside(ctl_points[-2:]) != begin_inside:
            bezier_path = concat([ctl_points_old[-2:], ctl_points])
            break

        ctl_points_old = ctl_points

    if bezier_path is None:
        raise ValueError("The path does not seem to intersect with the patch")

    bp = zip(bezier_path[::2], bezier_path[1::2])
    left, right = split_bezier_intersecting_with_closedpath(bp,
                                                            inside,
                                                            tolerence)
    if len(left) == 2:
        codes_left = [Path.LINETO]
        codes_right = [Path.MOVETO, Path.LINETO]
    elif len(left) == 3:
        codes_left = [Path.CURVE3, Path.CURVE3]
        codes_right = [Path.MOVETO, Path.CURVE3, Path.CURVE3]
    elif len(left) == 4:
        codes_left = [Path.CURVE4, Path.CURVE4, Path.CURVE4]
        codes_right = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]
    else:
        raise ValueError()

    verts_left = left[1:]
    verts_right = right[:]

    #i += 1

    if path.codes is None:
        path_in = Path(concat([path.vertices[:i], verts_left]))
        path_out = Path(concat([verts_right, path.vertices[i:]]))

    else:
        path_in = Path(concat([path.vertices[:iold], verts_left]),
                       concat([path.codes[:iold], codes_left]))

        path_out = Path(concat([verts_right, path.vertices[i:]]),
                        concat([codes_right, path.codes[i:]]))

    if reorder_inout and begin_inside == False:
        path_in, path_out = path_out, path_in

    return path_in, path_out





def inside_circle(cx, cy, r):
    r2 = r**2
    def _f(xy):
        x, y = xy
        return (x-cx)**2 + (y-cy)**2 < r2
    return _f



# quadratic bezier lines

def get_cos_sin(x0, y0, x1, y1):
    dx, dy = x1-x0, y1-y0
    d = (dx*dx + dy*dy)**.5
    return dx/d, dy/d

def check_if_parallel(dx1, dy1, dx2, dy2, tolerence=1.e-5):
    """ returns
       * 1 if two lines are parralel in same direction
       * -1 if two lines are parralel in opposite direction
       * 0 otherwise
    """
    theta1 = np.arctan2(dx1, dy1)
    theta2 = np.arctan2(dx2, dy2)
    dtheta = np.abs(theta1 - theta2)
    if  dtheta < tolerence:
        return 1
    elif np.abs(dtheta - np.pi) < tolerence:
        return -1
    else:
        return False


def get_parallels(bezier2, width):
    """
    Given the quadraitc bezier control points *bezier2*, returns
    control points of quadrativ bezier lines roughly parralel to given
    one separated by *width*.
    """

    # The parallel bezier lines constructed by following ways.
    #  c1 and  c2 are contol points representing the begin and end of the bezier line.
    #  cm is the middle point
    c1x, c1y = bezier2[0]
    cmx, cmy = bezier2[1]
    c2x, c2y = bezier2[2]

    parallel_test = check_if_parallel(c1x-cmx, c1y-cmy, cmx-c2x, cmy-c2y)

    if parallel_test == -1:
        warnings.warn("Lines do not intersect. A straight line is used instead.")
        #cmx, cmy = 0.5*(c1x+c2x), 0.5*(c1y+c2y)
        cos_t1, sin_t1 = get_cos_sin(c1x, c1y, c2x, c2y)
        cos_t2, sin_t2 = cos_t1, sin_t1
    else:
        # t1 and t2 is the anlge between c1 and cm, cm, c2.  They are
        # also a angle of the tangential line of the path at c1 and c2
        cos_t1, sin_t1 = get_cos_sin(c1x, c1y, cmx, cmy)
        cos_t2, sin_t2 = get_cos_sin(cmx, cmy, c2x, c2y)

    # find c1_left, c1_right which are located along the lines
    # throught c1 and perpendicular to the tangential lines of the
    # bezier path at a distance of width. Same thing for c2_left and
    # c2_right with respect to c2.
    c1x_left, c1y_left, c1x_right, c1y_right = \
              get_normal_points(c1x, c1y, cos_t1, sin_t1, width)
    c2x_left, c2y_left, c2x_right, c2y_right = \
              get_normal_points(c2x, c2y, cos_t2, sin_t2, width)

    # find cm_left which is the intersectng point of a line through
    # c1_left with angle t1 and a line throught c2_left with angle
    # t2. Same with cm_right.
    if parallel_test != 0:
        # a special case for a straight line, i.e., angle between two
        # lines are smaller than some (arbitrtay) value.
        cmx_left, cmy_left = \
                  0.5*(c1x_left+c2x_left), 0.5*(c1y_left+c2y_left)
        cmx_right, cmy_right = \
                   0.5*(c1x_right+c2x_right), 0.5*(c1y_right+c2y_right)
    else:
        cmx_left, cmy_left = \
                  get_intersection(c1x_left, c1y_left, cos_t1, sin_t1,
                                   c2x_left, c2y_left, cos_t2, sin_t2)

        cmx_right, cmy_right = \
                   get_intersection(c1x_right, c1y_right, cos_t1, sin_t1,
                                    c2x_right, c2y_right, cos_t2, sin_t2)

    # the parralel bezier lines are created with control points of
    # [c1_left, cm_left, c2_left] and [c1_right, cm_right, c2_right]
    path_left = [(c1x_left, c1y_left), (cmx_left, cmy_left), (c2x_left, c2y_left)]
    path_right = [(c1x_right, c1y_right), (cmx_right, cmy_right), (c2x_right, c2y_right)]

    return path_left, path_right



def make_wedged_bezier2(bezier2, length, shrink_factor=0.5):
    """
    Being similar to get_parallels, returns
    control points of two quadrativ bezier lines having a width roughly parralel to given
    one separated by *width*.
    """

    xx1, yy1 = bezier2[2]
    xx2, yy2 = bezier2[1]
    xx3, yy3 = bezier2[0]

    cx, cy = xx3, yy3
    x0, y0 = xx2, yy2

    dist = sqrt((x0-cx)**2 + (y0-cy)**2)
    cos_t, sin_t = (x0-cx)/dist, (y0-cy)/dist,

    x1, y1, x2, y2 = get_normal_points(cx, cy, cos_t, sin_t, length)

    xx12, yy12 = (xx1+xx2)/2., (yy1+yy2)/2.,
    xx23, yy23 = (xx2+xx3)/2., (yy2+yy3)/2.,

    dist = sqrt((xx12-xx23)**2 + (yy12-yy23)**2)
    cos_t, sin_t = (xx12-xx23)/dist, (yy12-yy23)/dist,

    xm1, ym1, xm2, ym2 = get_normal_points(xx2, yy2, cos_t, sin_t, length*shrink_factor)

    l_plus = [(x1, y1), (xm1, ym1), (xx1, yy1)]
    l_minus = [(x2, y2), (xm2, ym2), (xx1, yy1)]

    return l_plus, l_minus


def find_control_points(c1x, c1y, mmx, mmy, c2x, c2y):
    """ Find control points of the bezier line throught c1, mm, c2. We
    simply assume that c1, mm, c2 which have parameteric value 0, 0.5, and 1.
    """

    cmx = .5 * (4*mmx - (c1x + c2x))
    cmy = .5 * (4*mmy - (c1y + c2y))

    return [(c1x, c1y), (cmx, cmy), (c2x, c2y)]


def make_wedged_bezier2(bezier2, width, w1=1., wm=0.5, w2=0.):
    """
    Being similar to get_parallels, returns
    control points of two quadrativ bezier lines having a width roughly parralel to given
    one separated by *width*.
    """

    # c1, cm, c2
    c1x, c1y = bezier2[0]
    cmx, cmy = bezier2[1]
    c3x, c3y = bezier2[2]


    # t1 and t2 is the anlge between c1 and cm, cm, c3.
    # They are also a angle of the tangential line of the path at c1 and c3
    cos_t1, sin_t1 = get_cos_sin(c1x, c1y, cmx, cmy)
    cos_t2, sin_t2 = get_cos_sin(cmx, cmy, c3x, c3y)

    # find c1_left, c1_right which are located along the lines
    # throught c1 and perpendicular to the tangential lines of the
    # bezier path at a distance of width. Same thing for c3_left and
    # c3_right with respect to c3.
    c1x_left, c1y_left, c1x_right, c1y_right = \
              get_normal_points(c1x, c1y, cos_t1, sin_t1, width*w1)
    c3x_left, c3y_left, c3x_right, c3y_right = \
              get_normal_points(c3x, c3y, cos_t2, sin_t2, width*w2)




    # find c12, c23 and c123 which are middle points of c1-cm, cm-c3 and c12-c23
    c12x, c12y = (c1x+cmx)*.5, (c1y+cmy)*.5
    c23x, c23y = (cmx+c3x)*.5, (cmy+c3y)*.5
    c123x, c123y = (c12x+c23x)*.5, (c12y+c23y)*.5

    # tangential angle of c123 (angle between c12 and c23)
    cos_t123, sin_t123 = get_cos_sin(c12x, c12y, c23x, c23y)

    c123x_left, c123y_left, c123x_right, c123y_right = \
                get_normal_points(c123x, c123y, cos_t123, sin_t123, width*wm)


    path_left = find_control_points(c1x_left, c1y_left,
                                    c123x_left, c123y_left,
                                    c3x_left, c3y_left)
    path_right = find_control_points(c1x_right, c1y_right,
                                     c123x_right, c123y_right,
                                     c3x_right, c3y_right)

    return path_left, path_right




def make_path_regular(p):
    """
    fill in the codes if None.
    """
    c = p.codes
    if c is None:
        c = np.empty(p.vertices.shape[:1], "i")
        c.fill(Path.LINETO)
        c[0] = Path.MOVETO

        return Path(p.vertices, c)
    else:
        return p

def concatenate_paths(paths):
    """
    concatenate list of paths into a single path.
    """

    vertices = []
    codes = []
    for p in paths:
        p = make_path_regular(p)
        vertices.append(p.vertices)
        codes.append(p.codes)

    _path = Path(np.concatenate(vertices),
                 np.concatenate(codes))
    return _path



if 0:
    path = Path([(0, 0), (1, 0), (2, 2)],
                [Path.MOVETO, Path.CURVE3, Path.CURVE3])
    left, right = divide_path_inout(path, inside)
    clf()
    ax = gca()