## # # Most of the source code, except for the ones dealing with matplotlib are # a rough and direct translation from javascript code for bezier clipping. # This was written to accompany the ipython notebook viewable at # http://nbviewer.ipython.org/gist/hkrish/6527b9ad7cb7c2cc0cb0 # # The original javascript source code is available under MIT licence from # https://github.com/paperjs/paper.js/ # # The author disclaims copyright to this source code. # # /Harikrishnan Gopalakrishnan # /2015 # ## from math import * import matplotlib from matplotlib.path import Path import matplotlib.patches as patches import matplotlib.pyplot as plt import numpy as np from matplotlib.text import Text from operator import lt, gt class Point: x = 0 y = 0 def __init__(self, x, y): self.x = x self.y = y def __repr__(self): return str.format("< {}, {} >", self.x, self.y) def getSignedDistance(l1x, l1y, l2x, l2y, x, y): vx = l2x - l1x; vy = l2y - l1y; if (np.abs(vx) <= 1e-10): return l1x - x if vy >= 0 else x - l1x; m = vy / vx; # slope b = l1y - m * l1x; # y offset # Distance to the linear equation return (y - (m * x) - b) / np.sqrt(m * m + 1); def getConvexHull(dq0, dq1, dq2, dq3) : p0 = [ 0, dq0 ] p1 = [ 1. / 3, dq1 ] p2 = [ 2. / 3, dq2 ] p3 = [ 1, dq3 ] # Find signed distance of p1 and p2 from line [ p0, p3 ] dist1 = getSignedDistance(0, dq0, 1, dq3, 1. / 3, dq1) dist2 = getSignedDistance(0, dq0, 1, dq3, 2. / 3, dq2) flip = False hull = None; # Check if p1 and p2 are on the same side of the line [ p0, p3 ] if (dist1 * dist2 < 0) : # p1 and p2 lie on different sides of [ p0, p3 ]. The hull is a # quadrilateral and line [ p0, p3 ] is NOT part of the hull so we # are pretty much done here. # The top part includes p1, # we will reverse it later if that is not the case hull = [[p0, p1, p3], [p0, p2, p3]]; flip = dist1 < 0; else: # p1 and p2 lie on the same sides of [ p0, p3 ]. The hull can be # a triangle or a quadrilateral and line [ p0, p3 ] is part of the # hull. Check if the hull is a triangle or a quadrilateral. # Also, if at least one of the distances for p1 or p2, from line # [p0, p3] is zero then hull must at most have 3 vertices. pmax = None cross = 0 distZero = dist1 == 0 or dist2 == 0; if (abs(dist1) > abs(dist2)): pmax = p1; # apex is dq3 and the other apex point is dq0 vector dqapex -> # dqapex2 or base vector which is already part of the hull. cross = (dq3 - dq2 - (dq3 - dq0) / 3.) * (2 * (dq3 - dq2) - dq3 + dq1) / 3.; else: pmax = p2; # apex is dq0 in this case, and the other apex point is dq3 # vector dqapex -> dqapex2 or base vector which is already part # of the hull. cross = (dq1 - dq0 + (dq0 - dq3) / 3.) * (-2 * (dq0 - dq1) + dq0 - dq2) / 3.; # Compare cross products of these vectors to determine if the point # is in the triangle [ p3, pmax, p0 ], or if it is a quadrilateral. hull = [[p0, pmax, p3], [p0, p3]] if cross < 0 or distZero else [[p0, p1, p2, p3], [p0, p3]]; flip = dist1 < 0 if dist1 else dist2 < 0; if flip: hull.reverse(); return hull; def clipConvexHull(hullTop, hullBottom, dMin, dMax) : if (hullTop[0][1] < dMin): # Left of hull is below dMin, walk through the hull until it # enters the region between dMin and dMax return clipConvexHullPart(hullTop, True, dMin); elif (hullBottom[0][1] > dMax) : # Left of hull is above dMax, walk through the hull until it # enters the region between dMin and dMax return clipConvexHullPart(hullBottom, False, dMax); else : # Left of hull is between dMin and dMax, no clipping possible return hullTop[0][0]; def clipConvexHullPart(part, top, threshold) : px = part[0][0]; py = part[0][1]; for i in range(1, len(part)): qx = part[i][0]; qy = part[i][1]; if (qy >= threshold if top else qy <= threshold): return px + (threshold - py) * (qx - px) / (qy - py); px = qx; py = qy; return None; def addLocation(locations, t1, point1, t2, point2) : loc = [t1, point1, t2, point2]; # var loc = new CurveLocation(curve1, t1, point1, curve2, t2, point2); locations.append(loc); def subdivide(v, t = None) : p1x = v[0] p1y = v[1] c1x = v[2] c1y = v[3] c2x = v[4] c2y = v[5] p2x = v[6] p2y = v[7] if(t == None): t = 0.5; # Triangle computation, with loops unrolled. u = 1 - t # Interpolate from 4 to 3 points p3x = u * p1x + t * c1x p3y = u * p1y + t * c1y p4x = u * c1x + t * c2x p4y = u * c1y + t * c2y p5x = u * c2x + t * p2x p5y = u * c2y + t * p2y # Interpolate from 3 to 2 points p6x = u * p3x + t * p4x p6y = u * p3y + t * p4y p7x = u * p4x + t * p5x p7y = u * p4y + t * p5y # Interpolate from 2 points to 1 point p8x = u * p6x + t * p7x p8y = u * p6y + t * p7y; # We now have all the values we need to build the sub-curves: return [[p1x, p1y, p3x, p3y, p6x, p6y, p8x, p8y], [p8x, p8y, p7x, p7y, p5x, p5y, p2x, p2y]]; def getPart(v, frm, to) : if (frm > 0): v = subdivide(v, frm)[1]; # [1] right # Interpolate the parameter at 'to' in the new curve and cut there. if (to < 1): v = subdivide(v, (to - frm) / (1. - frm))[0]; # [0] left return v; def evaluate(v, t, typ) : EPSILON = 1e-10 TOLERANCE = 1e-5 p1x = v[0] p1y = v[1] c1x = v[2] c1y = v[3] c2x = v[4] c2y = v[5] p2x = v[6] p2y = v[7] tolerance = TOLERANCE x = 0. y = 0. # Handle special case at beginning / end of curve if (typ == 0 and (t < tolerance or t > 1 - tolerance)): isZero = t < tolerance; x = p1x if isZero else p2x; y = p1y if isZero else p2y; else : # Calculate the polynomial coefficients. cx = 3. * (c1x - p1x) bx = 3. * (c2x - c1x) - cx ax = p2x - p1x - cx - bx cy = 3. * (c1y - p1y) by = 3. * (c2y - c1y) - cy ay = p2y - p1y - cy - by if (typ == 0): # Calculate the curve point at parameter value t x = ((ax * t + bx) * t + cx) * t + p1x; y = ((ay * t + by) * t + cy) * t + p1y; else : # 1: tangent, 1st derivative # 2: normal, 1st derivative # 3: curvature, 1st derivative & 2nd derivative # Prevent tangents and normals of length 0: # http:#stackoverflow.com/questions/10506868/ if (t < tolerance and c1x == p1x and c1y == p1y or t > 1 - tolerance and c2x == p2x and c2y == p2y): x = c2x - c1x; y = c2y - c1y; elif (t < tolerance) : x = cx; y = cy; elif (t > 1 - tolerance) : x = 3. * (p2x - c2x); y = 3. * (p2y - c2y); else : # Simply use the derivation of the bezier function for both # the x and y coordinates: x = (3. * ax * t + 2 * bx) * t + cx; y = (3. * ay * t + 2 * by) * t + cy; pass if (typ == 3): # Calculate 2nd derivative, and curvature from there: # http:#cagd.cs.byu.edu/~557/text/ch2.pdf page#31 # k = |dx * d2y - dy * d2x| / (( dx^2 + dy^2 )^(3/2)) x2 = 6. * ax * t + 2. * bx; y2 = 6. * ay * t + 2. * by; return (x * y2 - y * x2) / pow(x * x + y * y, 3. / 2); # The normal is simply the rotated tangent: return Point(y, -x) if typ == 2 else Point(x, y); def addCurveIntersections(v1, v2, locations, curves=None, tMin=0., tMax=1., uMin=0., uMax=1., oldTDiff=1., reverse=False, recursion=0, recursionLimit=32, notes=[], data=None, tLimit=0.8) : # Avoid deeper recursion. # NOTE: @iconexperience determined that more than 20 recursions are # needed sometimes, depending on the tDiff threshold values further # below when determining which curve converges the least. He also # recommended a threshold of 0.5 instead of the initial 0.8 # See: https:#github.com/paperjs/paper.js/issues/565 if data <> None: prevSubdivide = data['subdivide'] if 'subdivide' in data else 0; if data['branchpoint']: prevSubdivide += 1 _instrProxy = {'reverse': reverse, 'step': recursion + 1, 'subdivide': prevSubdivide, 'limited': False, 'ixFound': False, 'relTDiff': 0, 'tDiff': 0, 'exit': False, 'children': [], 'branchpoint':False}; data['children'].append(_instrProxy) else: _instrProxy = None; if (recursion > recursionLimit): if data <> None: _instrProxy['limited'] = True; return; # Let P be the first curve and Q be the second TOLERANCE = 1e-5 q0x = v2[0] q0y = v2[1] q3x = v2[6] q3y = v2[7] tolerance = TOLERANCE epsilon = 1e-10 # Numerical.EPSILON, # Calculate the fat-line L for Q is the baseline l and two # offsets which completely encloses the curve P. d1 = getSignedDistance(q0x, q0y, q3x, q3y, v2[2], v2[3]) or 0 d2 = getSignedDistance(q0x, q0y, q3x, q3y, v2[4], v2[5]) or 0 factor = 3. / 4. if d1 * d2 > 0 else 4. / 9. dMin = factor * min(0, d1, d2) dMax = factor * max(0, d1, d2) # Calculate non-parametric bezier curve D(ti, di(t)) - di(t) is the # distance of P from the baseline l of the fat-line, ti is equally # spaced in [0, 1] dp0 = getSignedDistance(q0x, q0y, q3x, q3y, v1[0], v1[1]) dp1 = getSignedDistance(q0x, q0y, q3x, q3y, v1[2], v1[3]) dp2 = getSignedDistance(q0x, q0y, q3x, q3y, v1[4], v1[5]) dp3 = getSignedDistance(q0x, q0y, q3x, q3y, v1[6], v1[7]) tMinNew = 0. tMaxNew = 0. tDiff = 0. # NOTE: the recursion threshold of 4 is needed to prevent issue #571 # from occurring: https:#github.com/paperjs/paper.js/issues/571 if (q0x == q3x and uMax - uMin <= epsilon and recursion > 4): # The fatline of Q has converged to a point, the clipping is not # reliable. Return the value we have even though we will miss the # precision. tMaxNew = tMinNew = (tMax + tMin) / 2.; tDiff = 0; if data <> None: _instrProxy['tDiff'] = 0; _instrProxy['relTDiff'] = 0; else : # Get the top and bottom parts of the convex-hull hull = getConvexHull(dp0, dp1, dp2, dp3) top = hull[0] bottom = hull[1] tMinClip = None tMaxClip = None # Clip the convex-hull with dMin and dMax tMinClip = clipConvexHull(top, bottom, dMin, dMax); top.reverse(); bottom.reverse(); tMaxClip = clipConvexHull(top, bottom, dMin, dMax); # No intersections if one of the tvalues are null or 'undefined' if (tMinClip == None or tMaxClip == None): if data <> None: _instrProxy['exit'] = True; return; # Clip P with the fatline for Q v1d = v1 v1 = getPart(v1, tMinClip, tMaxClip); tDiff = tMaxClip - tMinClip; if curves <> None: curves.append({'before': v1d, 'after': v1, 'clipper': v2, 'notes':notes, 'tdiff': tDiff, 'ix': 0, 'subdivide': None}) # tMin and tMax are within the range (0, 1). We need to project it # to the original parameter range for v2. tMinNew = tMax * tMinClip + tMin * (1 - tMinClip); tMaxNew = tMax * tMaxClip + tMin * (1 - tMaxClip); if data <> None: _instrProxy['tDiff'] = tMaxNew - tMinNew; _instrProxy['relTDiff'] = tDiff; pass # Check if we need to subdivide the curves if (oldTDiff > tLimit and tDiff > tLimit): # Subdivide the curve which has converged the least. if data <> None: _instrProxy['branchpoint'] = True if (tMaxNew - tMinNew > uMax - uMin): parts = subdivide(v1, 0.5) if curves <> None: curves[-1]['subdivide'] = [parts[1][0], parts[1][1]] t = tMinNew + (tMaxNew - tMinNew) / 2. recursion += 1 addCurveIntersections(v2, parts[0], locations, curves, uMin, uMax, tMinNew, t, tDiff, not reverse, recursion, recursionLimit, notes=['left'], data=_instrProxy, tLimit=tLimit); addCurveIntersections(v2, parts[1], locations, curves, uMin, uMax, t, tMaxNew, tDiff, not reverse, recursion, recursionLimit, notes=['right'], data=_instrProxy, tLimit=tLimit); else : parts = subdivide(v2, 0.5) if curves <> None: curves[-1]['subdivide'] = [parts[1][0], parts[1][1]] t = uMin + (uMax - uMin) / 2.; recursion += 1 addCurveIntersections(parts[0], v1, locations, curves, uMin, t, tMinNew, tMaxNew, tDiff, not reverse, recursion, recursionLimit, notes=['left'], data=_instrProxy, tLimit=tLimit); addCurveIntersections(parts[1], v1, locations, curves, t, uMax, tMinNew, tMaxNew, tDiff, not reverse, recursion, recursionLimit, notes=['right'], data=_instrProxy, tLimit=tLimit); pass pass elif (max(uMax - uMin, tMaxNew - tMinNew) < tolerance) : # We have isolated the intersection with sufficient precision t1 = tMinNew + (tMaxNew - tMinNew) / 2. t2 = uMin + (uMax - uMin) / 2. if data <> None: _instrProxy['ixFound'] = True if curves <> None: curves[-1]['ix'] = len(locations) + 1; if (reverse): addLocation(locations, t2, evaluate(v2, t2, 0), t1, evaluate(v1, t1, 0)); else : addLocation(locations, t1, evaluate(v1, t1, 0), t2, evaluate(v2, t2, 0)); pass pass else: recursion += 1 addCurveIntersections(v2, v1, locations, curves, uMin, uMax, tMinNew, tMaxNew, tDiff, not reverse, recursion, recursionLimit, notes, _instrProxy, tLimit=tLimit); pass def getFatline(v): "This method returns a tuple containing maximum and minimum offset width for the 'fatline'." # Starting point of the curve q0x = v[0] q0y = v[1] # End point of the curve q3x = v[6] q3y = v[7] # Calculate the fat-line L, for Q is the baseline l and two # offsets which completely encloses the curve P. d1 = getSignedDistance(q0x, q0y, q3x, q3y, v[2], v[3]) or 0 d2 = getSignedDistance(q0x, q0y, q3x, q3y, v[4], v[5]) or 0 factor = 3. / 4. if d1 * d2 > 0 else 4. / 9. dMin = factor * min(0, d1, d2) dMax = factor * max(0, d1, d2) # The width of the 'fatline' is |dMin| + |dMax| return (dMin, dMax) def getCrossings(v1, v2, steps=32, fig=None, subplots=None, size=None, labels=None): loc = [] curves = [] addCurveIntersections(v1, v2, locations=loc, curves=curves, recursionLimit=steps) colors = ['#e74c3c', '#34495e'] labels = ['$Curve\ A$', '$Curve\ B$'] if labels == None else labels; ty = 0 tx = 0 boldText = matplotlib.font_manager.FontProperties(weight='bold') if subplots <> None: count = min(len(subplots), len(curves)) for i in range(count): clip = curves[i] ax = subplots[i] if size == None: size = drawCurve(ax, v1, setsize=True, color='#eeeeee')[0] else: ax.set_xlim(left=-size-10, right=size+10) ax.set_ylim(bottom=-size-10, top=size+10) drawCurve(ax, v1, color='#eeeeee') drawCurve(ax, v2, color='#eeeeee') drawCurve(ax, clip['before'], color='#aaaaaa') drawCurve(ax, clip['after'], color=colors[(i+1)%2], label=labels[(i+1)%2]) drawCurve(ax, clip['clipper'], color=colors[i%2], fatdist=getFatline(clip['clipper']), label=labels[i%2]) tx = -size-5 ty = size ix = clip['ix'] if ix > 0: ax.text(tx, ty, str.format('Intersection {} found', ix), color=colors[0], fontproperties=boldText) ty -= 9 notes = clip['notes'] ax.text(tx, ty, '$Length\ after\ clip$: ' + str.format('${:3.0f}\%\ of\ prev.$', clip['tdiff']*100)) if len(notes) > 0: ty -= 9 ax.text(tx, ty, '$Side$: ' + str(notes[0])) sub = clip['subdivide'] if sub <> None: ty -= 18 ax.text(tx, ty, str.format('$Clipped\ <20\%\ of\ the\ length.$\n$Subdividing...$', ix), color=colors[0]) ax.plot([sub[0]], [sub[1]], 'x', markeredgecolor='w', markeredgewidth=3, markersize=7.0) ax.plot([sub[0]], [sub[1]], 'x', markeredgecolor='#3498db', markeredgewidth=2, markersize=6.0) for i in range(count, len(subplots)): fig.delaxes(subplots[i]) return loc def accumulateTDiffs(dat, tDiffs): "Just collect the data we need at the moment." tDiffs.append((dat['step'], dat['subdivide'], dat['branchpoint'], dat['exit'], dat['ixFound'], dat['tDiff'],dat['relTDiff'])) children = dat['children'] dat.pop('children') dat['name'] = str.format("<{}, {}>", dat['step'], dat['subdivide']) if len(children) == 0: return Tree(node=dat, children=[]) chArray = [] for f in children: chArray.append(accumulateTDiffs(f, tDiffs)) return Tree(node=dat, children=chArray) def getCrossingsInstr(v1, v2, tLimit=0.8): loc = [] data = {'reverse': False, 'step': 0, 'subdivide': 0, 'limited': False, 'ixFound': False, 'relTDiff': 1, 'tDiff': 1, 'exit': False, 'children': [], 'branchpoint':False} addCurveIntersections(v1, v2, locations=loc, data=data, tLimit=tLimit) tDiffs = [] tree = buchheim(accumulateTDiffs(data, tDiffs)) return (loc, tDiffs, tree) def getCrossingsPlain(v1, v2, tLimit=0.8): loc = [] addCurveIntersections(v1, v2, locations=loc, tLimit=tLimit) return loc # Helper methods def drawCurve(ax, v, color="#000000", fatcolor=None, lw=2, fatdist=(0,0), setsize=False, label='', annotate=False, points=True): v = np.array(v) codes = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4] xy = v.reshape((2,4), order='F') verts = xy.T path = Path(verts, codes) ret = [0, 0] # plot if setsize: x, y = xy size = np.max((np.max(x) - np.min(x), np.max(y) - np.min(y))) / 2. ax.set_xlim(left=-size-10, right=size+10) ax.set_ylim(bottom=-size-10, top=size+10) ret[0] = size if color == None: return ax if fatdist != (0, 0): colorconv = matplotlib.colors.ColorConverter() fatcolor = color if fatcolor == None else fatcolor; ret[1] = drawFatLine(ax, v, fatdist, fatcolor, fatcolor) patch = patches.PathPatch(path, facecolor='none', edgecolor=color, lw=lw) ax.add_patch(patch) if points: x, y = zip(*verts) ax.plot(x, y, ls='--', marker='o', markeredgecolor='w', color=color, label=label) if annotate: for i, v in enumerate(verts): ax.text(v[0], v[1], str.format('$V_{}$', i), size=15) return ret def drawLine(ax, v, color="#000000", lw=2, setsize=False, zorder=1): codes = [Path.MOVETO, Path.LINETO] xy = v.reshape((2,2), order='F') verts = xy.T path = Path(verts, codes) # plot if setsize: x, y = xy size = np.max((np.max(x) - np.min(x), np.max(y) - np.min(y))) / 2. ax.set_xlim(left=-size-10, right=size+10) ax.set_ylim(bottom=-size-10, top=size+10) patch = patches.PathPatch(path, facecolor='none', edgecolor=color, lw=lw, zorder=zorder) ax.add_patch(patch) return ax def drawPoly(ax, v, ecolor='none', fcolor='#000000', lw=1, setsize=False): ncodes = (len(v) / 2) - 1 codes = [Path.MOVETO].extend( [Path.LINETO for i in range(ncodes)]); xy = v.reshape((2,-1), order='F') verts = xy.T path = Path(verts, codes) # plot if setsize: x, y = xy size = np.max((np.max(x) - np.min(x), np.max(y) - np.min(y))) / 2. ax.set_xlim(left=-size-10, right=size+10) ax.set_ylim(bottom=-size-10, top=size+10) patch = patches.PathPatch(path, facecolor=fcolor, edgecolor=ecolor, lw=lw, alpha=0.3) ax.add_patch(patch) return ax def drawRect(ax, x, y, width, height, ecolor='none', fcolor='#000000', lw=1): patch = patches.Rectangle(xy=(x, y), width=width, height=height, facecolor=fcolor, edgecolor=ecolor, lw=lw) ax.add_patch(patch) return ax def drawFatLine(ax, crv, fatDist, ecolor='#000000', fcolor=(0,0,0,0.3)): baseline = np.concatenate([crv[0:2], crv[6:8]]) drawLine(ax, baseline, ecolor, lw=0.5) baselineSlope = np.arctan2(baseline[3] - baseline[1], baseline[2] - baseline[0]) signum = 1 baselineSlope = (baselineSlope + 2*np.pi) % (2*np.pi) if baselineSlope <= 3*np.pi/2. and baselineSlope > np.pi/2.: signum *= -1 cosa = np.cos(baselineSlope) sina = -np.sin(baselineSlope) fatline = np.array([]) fatDist = np.add(fatDist, (-0.5, 0.5)) # Offset the line thickness lines = [] for d in fatDist: ln = baseline.copy() mat = np.matrix([[cosa, sina, 0.], [-sina, cosa, 0.], [0., 0., 1.]]) coeffs = np.array((mat * np.matrix([[0], [signum * d], [1]]))[0:2, :]).reshape(-1) mat = np.identity(3, dtype=np.float64) mat[0][2] = coeffs[0] mat[1][2] = coeffs[1] m = np.matrix(np.vstack((ln.reshape((2,2), order='F'), np.ones(2)))) ln = np.array((mat * m)[0:2,:]).reshape((-1), order='F') drawLine(ax, ln, ecolor, lw=0.5) fatline = np.hstack((fatline, ln)) lines.append(ln) fatlineVert = fatline.reshape((2, 4), order='F').T fatlineVert[[2,3], :] = fatlineVert[[3,2], :] drawPoly(ax, fatlineVert.reshape((-1)), fcolor=fcolor) return lines def rotate(v, angle): cosa = np.cos(angle) sina = np.sin(angle) mat = np.matrix([[cosa, sina, 0.], [-sina, cosa, 0.], [0., 0., 1.]]) ncols = len(v) / 2 verts = np.matrix(np.vstack((v.reshape((2,-1), order='F'), np.ones(ncols)))) verts = np.array((mat * verts)[0:2, :]).reshape(-1, order='F') return verts # Below are some methods, which just visualises various information def drawNonParamCurveI(fig, gs, curve1r, curve2, fatDist, nonParamB, chull, tclips): ax1 = plt.subplot(gs[0, 0]) ax2 = plt.subplot(gs[0, 1]) titles = ['non-parametric version of $B$'] ax1.set_xlabel('px') ax1.set_ylabel('px') ax1.set_title('$B$ in cartesian space') ax2.set_xlabel('$t_i$') ax2.set_ylabel('distance from $L$') ax2.set_title('Non-parametric version of $B$') # plot the curves colorconv = matplotlib.colors.ColorConverter() drawFatLine(ax1, curve1r, fatDist=fatDist, fcolor=colorconv.to_rgba('#e74c3c', alpha=0.3)) drawCurve(ax1, curve2, '#34495e', setsize=True, annotate=True) sp = curve1r[0:2]; ep = curve1r[6:]; txy=sp + 3. * (ep - sp) / 4. ax1.text(txy[0], txy[1]+2.5, "$L$", ha="center", color='black', va="center",size=20) ax1.text(txy[0], txy[1]+fatDist[0]-3.5, "$d_{min}$", ha="center", color='black', va="center",size=20) ax1.text(txy[0], txy[1]+fatDist[1]+3.5, "$d_{max}$", ha="center", color='black', va="center",size=20) # plot non-parametric curve drawPoly(ax2, chull, ecolor='#3498db', fcolor=colorconv.to_rgba('#3498db', alpha=0.3)) drawCurve(ax2, nonParamB, '#34495e', annotate=True) ax2.axhline(0, color='k') ax2.axhline(fatDist[0], color='#e74c3c') ax2.text(0.1, fatDist[0], "$d_{min}$", ha="center", color='black', size=20) ax2.axhline(fatDist[1], color='#e74c3c') ax2.text(0.1, fatDist[1], "$d_{max}$", ha="center", color='black', size=20) # plot tMin and tMax nonParamBxy = nonParamB.reshape((2,-1), order='F') ymin = min(nonParamBxy[1]) ymax = max(nonParamBxy[1]) tMin, tMax = tclips tMinL = np.array([tMin, ymin, tMin, fatDist[0]]) drawLine(ax2, tMinL, color='#2ecc71') tMaxL = np.array([tMax, ymin, tMax, fatDist[1]]) drawLine(ax2, tMaxL, color='#2ecc71') dx = ymin + (ymax - ymin)/20. ax2.text(tMin, dx, "$t_{min}$", ha="center", color='black', size=20) ax2.text(tMax, dx, "$t_{max}$", ha="center", color='black', size=20) ax2.plot([tMin], [fatDist[0]], 'o', markeredgecolor='#2ecc71', color='#e74c3c') ax2.plot([tMax], [fatDist[1]], 'o', markeredgecolor='#2ecc71', color='#e74c3c') ax2.set_xlim(left=0., right=1.) ax2.set_ylim(bottom=ymin, top=ymax) # plot the same points on curve B also curve2P = getPart(curve2, tMin, tMax) drawCurve(ax1, curve2P, '#3498db', points=False) pMin = evaluate(curve2, tMin, 0) ax1.plot([pMin.x], [pMin.y], 'o', markeredgecolor='#2ecc71', color='#e74c3c') ax1.text(pMin.x + 1, pMin.y, "$t_{min}$", ha="left", color='black', size=20) pMax = evaluate(curve2, tMax, 0) ax1.plot([pMax.x], [pMax.y], 'o', markeredgecolor='#2ecc71', color='#e74c3c') ax1.text(pMax.x + 1, pMax.y, "$t_{max}$", ha="left", color='black', size=20) hull_patch = patches.Patch(color='#3498db', label='convex-hull') plt.legend(handles=[hull_patch], loc=0) plt.tight_layout(); #----------------------------------------------------------------------------- # Test data generation def genBezTestCasesBatch(nCases, slope=(0.,np.pi/2), size=(50., 100.), spread=(5., 100.)): """Generates a number of random bezier curves that roughly fits a centerline (slope) with defined size (pixels) and spread (distance of points of the curve spread from the centerline)""" # Generate a set of random values for x and y in the range of [-0.5, 0.5) randNums = np.random.rand(2 * nCases, 8) - 0.5 xvals, yvals = (randNums[:,:4], randNums[:,4:]) del(randNums) # Arrange the random values as x, y pairs for each curve # xvals.sort(axis=1) # xvals = np.hstack([np.repeat(-0.45, 2 * nCases).reshape((-1,1)), # xvals, # np.repeat(0.45, 2 * nCases).reshape((-1,1))]) # randNums = np.random.rand(2 * nCases, 6) - 0.5 # xvals, yvals = (randNums[:,:2], randNums[:,2:]) # del(randNums) # # Arrange the random values as x, y pairs for each curve # xvals.sort(axis=1) # xvals = np.hstack([np.repeat(-0.45, 2 * nCases).reshape((-1,1)), # xvals, # np.repeat(0.45, 2 * nCases).reshape((-1,1))]) verts = np.hstack([xvals[:,np.newaxis], yvals[:,np.newaxis], np.ones(xvals.shape)[:,np.newaxis]]) verts = verts.transpose((0,2,1)) # Initialise the transformation matrix transMat = np.ones(18 * nCases).reshape((-1,3,3)) arr = np.array([(i, j) for i in range(3) for j in range(3) if i <> j]) arr = np.vsplit(arr.T,1)[0] transMat[:, arr[0], arr[1]] = 0. del(arr) # Transformation: rotation slopeMat = transMat.copy() angle = slope[0] + np.random.random(2 * nCases) * (slope[1] - slope[0]) cosa = np.cos(angle) sina = np.sin(angle) slopeMat[:, 0, 0] = slopeMat[:, 1, 1] = cosa slopeMat[:, 0, 1] = sina slopeMat[:, 1, 0] = -sina # Transformation: Scale sizes = size[0] + np.random.random(2 * nCases) * (size[1] - size[0]) spreads = spread[0] + np.random.random(2 * nCases) * (spread[1] - spread[0]) transMat[:, 0, 0] = sizes transMat[:, 1, 1] = spreads transVerts = [] # Apply the transformation and return for i, vert in enumerate(verts): m = np.dot(transMat[i], slopeMat[i]) a = np.dot(vert, m) transVerts.append(a[:,:2].reshape(-1)) transVerts = np.array(transVerts).reshape(-1, 2, 8) return transVerts def genBezTestCases(nCases, slope=(0.,np.pi/2), size=(50., 100.), spread=(5., 100.), ix=(1, 9)): """This method makes sure we have enough number of intersections in our test cases, and discards the ones that does not fit the criteria""" count = 0 cases = [] ixMin, ixMax = ix if ixMax < ixMin: ixMax, ixMin = (ixMin, ixMax) # if ixMin > 4: print 'Limiting the intersections to 4! It will take a long time to generate.' # ixMin = 4; batchLen = nCases if nCases > 1000 else 1000 while count < nCases: batch = genBezTestCasesBatch(batchLen, slope, size, spread) for b in batch: ixc = len(getCrossingsPlain(b[0], b[1])) if ixc >= ixMin and ixc <= ixMax: cases.append(b) count += 1 # print str(count) + ", ", if count >= nCases: break; return cases def drawClippingTree(fig, ax, tree, color='#3498db', ecolor=None, r=10, figsize=6, zorder=1, yshift=0, labels=None): "Visualises how the clipping algorithm recurves and subdivides using a tree." if labels == None: labels = {} labels['o'] = None labels['>'] = None labels['x'] = None labels['s'] = None wid, hei, span, depth = drawTree(ax, tree, r=r, color=color, ecolor=ecolor, zorder=1, yshift=yshift, labels=labels) drawTreeEdges(ax, tree, r=r, color=color, ecolor=ecolor, zorder=1, yshift=yshift) hei += yshift ax.set_xlim(left=-1, right=wid+1) ax.set_ylim(top=hei+1, bottom=-1) # resize the plot maxDepth = depth maxSpan = span maxD, minD = (max(span, depth), min(span, depth)) aspect = minD / (1. * maxD) maxD = figsize * 1. minD = aspect * maxD span, depth = (maxD, minD) if span > depth else (minD, maxD) fig.set_size_inches((depth, span)) return (wid, hei-yshift, maxSpan, maxDepth) def drawStackedTDiffs(fig, ax, tDiffs, maxY, minY=0.01): spacing = 1. maxDepth = max(tDiffs)[0] tDiffForDepth = [[] for x in range(maxDepth+1)] # Palette generated with the excellent seaborn package # http://stanford.edu/~mwaskom/software/seaborn/index.html palette = ["#66c2a5", "#fa8d62", "#8d9fca", "#e68ac3", "#a6d853", "#fed82f", "#e4c394", "#b3b3b3", "#66c2a5", "#fa8d62", "#8d9fca", "#e68ac3", "#a6d853", "#fed82f", "#e4c394", "#b3b3b3", "#66c2a5", "#fa8d62", "#8d9fca", "#e68ac3", "#a6d853", "#fed82f", "#e4c394", "#b3b3b3", "#66c2a5", "#fa8d62", "#8d9fca", "#e68ac3", "#a6d853", "#fed82f"] paletteDesat = ["#86a199", "#c5a497", "#a2a8b5", "#c6aabb", "#9aa982", "#b6aa77", "#c8beb0", "#b3b3b3", "#86a199", "#c5a497", "#a2a8b5", "#c6aabb", "#9aa982", "#b6aa77", "#c8beb0", "#b3b3b3", "#86a199", "#c5a497", "#a2a8b5", "#c6aabb", "#9aa982", "#b6aa77", "#c8beb0", "#b3b3b3", "#86a199", "#c5a497", "#a2a8b5", "#c6aabb", "#9aa982", "#b6aa77"] # np.random.shuffle(palette) minYaxes = 0 maxY = maxY if maxY > minYaxes else minYaxes + 1 ax.set_ylim(top=maxY, bottom=0) for tdiff in tDiffs: dep = tdiff[0] tDiffForDepth[dep].append(tdiff) sumTDiffs = map(lambda x: reduce(lambda tsum, y: tsum+y[5], x, 0), tDiffForDepth) maxtDiff = max(sumTDiffs) yPerT = 1 #/ (maxtDiff * 1.) blocks = [0] lastC = 0 for dep, tdep in enumerate(tDiffForDepth): totalT = sumTDiffs[dep] ybot = 0 diffB = [] # Plot each stack while keeping a tab on branching and # cancelled recursion as a diff for j, t in enumerate(tdep): bHei = minY + t[5] * yPerT color = palette[blocks[j]] if not t[3] else 'k' color = color if not t[2] else paletteDesat[blocks[j]] drawRect(ax, dep*spacing-spacing/2, ybot, spacing, bHei, fcolor=color, ecolor='k', lw=0.5) ax.plot(dep*spacing-spacing/2, ybot+bHei, marker='o', markersize=0) ybot+= bHei if t[2]: diffB.append(('+', j)) elif t[3]: diffB.append(('-', j)) pass # print 'diff ', diffB adjB = 0 for d in diffB: if d[0] == '+': idx = d[1] + adjB blocks[idx:idx+1] = [lastC+1, lastC+2] lastC = (lastC + 2) % (len(palette)-2) adjB += 1 elif d[0] == '-': idx = d[1] + adjB blocks[idx:idx+1] = [] adjB -= 1 pass ##============================================================================== # The following code draws a tree in O(n) time. # From http://billmill.org/pymag-trees/ # https://github.com/llimllib/pymag-trees/ class Tree: def __init__(self, node={}, children=[]): self.node = node self.width = len(node['name']) self.children = children def __str__(self): return "%s" % (self.node) def __repr__(self): return "%s" % (self.node) def __getitem__(self, key): if isinstance(key, int) or isinstance(key, slice): return self.children[key] if isinstance(key, str): for child in self.children: if child.node == key: return child def __iter__(self): return self.children.__iter__() def __len__(self): return len(self.children) class DrawTree(object): def __init__(self, tree, parent=None, depth=0, number=1): self.x = -1. self.y = depth self.tree = tree self.children = [DrawTree(c, self, depth+1, i+1) for i, c in enumerate(tree.children)] self.parent = parent self.thread = None self.mod = 0 self.ancestor = self self.change = self.shift = 0 self._lmost_sibling = None #this is the number of the node in its group of siblings 1..n self.number = number def left(self): return self.thread or len(self.children) and self.children[0] def right(self): return self.thread or len(self.children) and self.children[-1] def lbrother(self): n = None if self.parent: for node in self.parent.children: if node == self: return n else: n = node return n def get_lmost_sibling(self): if not self._lmost_sibling and self.parent and self != \ self.parent.children[0]: self._lmost_sibling = self.parent.children[0] return self._lmost_sibling lmost_sibling = property(get_lmost_sibling) def __str__(self): return "%s: x=%s mod=%s" % (self.tree, self.x, self.mod) def __repr__(self): return self.__str__() def buchheim(tree): dt = firstwalk(DrawTree(tree)) min = second_walk(dt) if min < 0: third_walk(dt, -min) return dt def third_walk(tree, n): tree.x += n for c in tree.children: third_walk(c, n) def firstwalk(v, distance=1.): if len(v.children) == 0: if v.lmost_sibling: v.x = v.lbrother().x + distance else: v.x = 0. else: default_ancestor = v.children[0] for w in v.children: firstwalk(w) default_ancestor = apportion(w, default_ancestor, distance) # print "finished v =", v.tree, "children" execute_shifts(v) midpoint = (v.children[0].x + v.children[-1].x) / 2 ell = v.children[0] arr = v.children[-1] w = v.lbrother() if w: v.x = w.x + distance v.mod = v.x - midpoint else: v.x = midpoint return v def apportion(v, default_ancestor, distance): w = v.lbrother() if w is not None: #in buchheim notation: #i == inner; o == outer; r == right; l == left; r = +; l = - vir = vor = v vil = w vol = v.lmost_sibling sir = sor = v.mod sil = vil.mod sol = vol.mod while vil.right() and vir.left(): vil = vil.right() vir = vir.left() vol = vol.left() vor = vor.right() vor.ancestor = v shift = (vil.x + sil) - (vir.x + sir) + distance if shift > 0: move_subtree(ancestor(vil, v, default_ancestor), v, shift) sir = sir + shift sor = sor + shift sil += vil.mod sir += vir.mod sol += vol.mod sor += vor.mod if vil.right() and not vor.right(): vor.thread = vil.right() vor.mod += sil - sor else: if vir.left() and not vol.left(): vol.thread = vir.left() vol.mod += sir - sol default_ancestor = v return default_ancestor def move_subtree(wl, wr, shift): subtrees = wr.number - wl.number # print wl.tree, "is conflicted with", wr.tree, 'moving', subtrees, 'shift', shift #print wl, wr, wr.number, wl.number, shift, subtrees, shift/subtrees wr.change -= shift / subtrees wr.shift += shift wl.change += shift / subtrees wr.x += shift wr.mod += shift def execute_shifts(v): shift = change = 0 for w in v.children[::-1]: # print "shift:", w, shift, w.change w.x += shift w.mod += shift change += w.change shift += w.shift + change def ancestor(vil, v, default_ancestor): #the relevant text is at the bottom of page 7 of #"Improving Walker's Algorithm to Run in Linear Time" by Buchheim et al, (2002) #http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.16.8757&rep=rep1&type=pdf if vil.ancestor in v.parent.children: return vil.ancestor else: return default_ancestor def second_walk(v, m=0, depth=0, min=None): v.x += m v.y = depth if min is None or v.x < min: min = v.x for w in v.children: min = second_walk(w, m + v.mod, depth+1, min) return min def drawTree(ax, root, depth=0, color='b', ecolor=None, r=10, zorder=1, yshift=0, labels={'o': None, '>': None, 'x': None, 's': None}): """Draw the nodes of the *clipping algorithm progression as a tree*""" rh = 1. rw = 1. ew = 1.5 siz = r shape = 'o' ecolor = 'w' if ecolor == None else ecolor node = root.tree.node label1 = '$t_{diff} \leq t_{limit}\ ({\\bf good\ clip})$' label2 = '' if node['branchpoint']: shape='>' label1 = '$t_{diff} > t_{limit}\ ({\\bf subdivide})$' elif node['exit']: shape='x' label1 = '' label2 = '$degenerate ({\\bf no\ intersections})$' ew = 3 elif node['ixFound']: shape='s' label1 = '$\\bf t_{diff} < \epsilon\ (isolated\ an\ intersection)$' else: siz = 2 * r / 3. if labels[shape] <> None: label1 = label2 = '' plt1, = ax.plot(depth * rh, root.x * rw + yshift, marker=shape, color=color, markersize=siz, markerfacecolor=color, markeredgecolor=ecolor, markeredgewidth=ew, zorder=zorder, label=label1) if shape == 'x': plt2, = ax.plot(depth * rh, root.x * rw + yshift, marker=shape, markersize=r-1, markerfacecolor=color, color=color, markeredgecolor=color, markeredgewidth=ew-1.5, zorder=zorder, label=label2) # Save one of each kind of node, we will use this later to create a legend. if shape=='o' and labels['o'] == None: labels['o'] = plt1 elif shape=='>' and labels['>'] == None: labels['>'] = plt1 elif shape=='x' and labels['x'] == None: labels['x'] = plt2 elif shape=='s' and labels['s'] == None: labels['s'] = plt1 if len(root.children) == 0: return (depth * rh, root.x * rw, 1, depth); maxWid = 0 maxHei = 0 maxDepth = 0 maxSpan = 0 for child in root.children: wid, hei, span, dep = drawTree(ax, child, depth+1, color=color, ecolor=ecolor, r=r, zorder=zorder, yshift=yshift, labels=labels) if hei > maxHei: maxHei = hei if wid > maxWid: maxWid = wid if dep > maxDepth: maxDepth = dep maxSpan += span return (maxWid, maxHei, maxSpan, maxDepth); def drawTreeEdges(ax, root, depth=0, color='k', ecolor=None, r=10, zorder=1, yshift=0): """Draw the edges of the *clipping algorithm progression as a tree*""" rh = 1. rw = 1. ew = 1.5 for child in root.children: lineVerts = np.array([depth * rh, root.x * rw + yshift, (depth+1) * rh, child.x * rw + yshift]) if ecolor <>None: drawLine(ax, lineVerts, color=ecolor, lw=2 * ew, zorder=zorder) drawLine(ax, lineVerts, color=color, lw=1, zorder=zorder) drawTreeEdges(ax, child, depth+1, color=color, ecolor=ecolor, r=r, zorder=zorder, yshift=yshift) pass; def getSpanDepth(root, depth=0): if len(root.children) == 0: return (1, depth); maxDepth = 0 maxSpan = 0 for child in root.children: span, dep = getSpanDepth(child, depth+1) return maxSpan, maxDepth # Display some values as HTML table in the IPython notebook. class ListTable(list): """ Overridden list class which takes a 2-dimensional list of the form [[1,2,3],[4,5,6]], and renders an HTML Table in IPython Notebook. """ def _repr_html_(self): html = [""] for row in self: html.append("