1# This file is part of the Astrometry.net suite. 2# Licensed under a 3-clause BSD style license - see LICENSE 3from math import pi 4import numpy as np 5 6def estimate_mode(img, lo=None, hi=None, plo=1, phi=70, bins1=30, 7 flo=0.5, fhi=0.8, bins2=30, 8 return_fit=False, raiseOnWarn=False): 9 # Estimate sky level by: compute the histogram within [lo,hi], fit 10 # a parabola to the log-counts, return the argmax of that parabola. 11 12 # Coarse bin to find the peak (mode) 13 if lo is None: 14 lo = np.percentile(img, plo) 15 if hi is None: 16 hi = np.percentile(img, phi) 17 18 binedges1 = np.linspace(lo, hi, bins1+1) 19 counts1,e = np.histogram(img.ravel(), bins=binedges1) 20 bincenters1 = binedges1[:-1] + (binedges1[1]-binedges1[0])/2. 21 maxbin = np.argmax(counts1) 22 maxcount = counts1[maxbin] 23 mode = bincenters1[maxbin] 24 25 # Search for bin containing < {flo,fhi} * maxcount 26 ilo = maxbin 27 while ilo > 0: 28 ilo -= 1 29 if counts1[ilo] < flo*maxcount: 30 break 31 ihi = maxbin 32 while ihi < bins1-1: 33 ihi += 1 34 if counts1[ihi] < fhi*maxcount: 35 break 36 37 lo = bincenters1[ilo] 38 hi = bincenters1[ihi] 39 40 binedges = np.linspace(lo, hi, bins2) 41 counts,e = np.histogram(img.ravel(), bins=binedges) 42 bincenters = binedges[:-1] + (binedges[1]-binedges[0])/2. 43 44 b = np.log10(np.maximum(1, counts)) 45 46 xscale = 0.5 * (hi - lo) 47 x0 = (hi + lo) / 2. 48 x = (bincenters - x0) / xscale 49 50 A = np.zeros((len(x), 3)) 51 A[:,0] = 1. 52 A[:,1] = x 53 A[:,2] = x**2 54 res = np.linalg.lstsq(A, b, rcond=None) 55 X = res[0] 56 mx = -X[1] / (2. * X[2]) 57 mx = (mx * xscale) + x0 58 59 warn = None 60 61 if not (mx > lo and mx < hi): 62 if raiseOnWarn: 63 raise ValueError('sky estimate not bracketed by peak: lo %f, sky %f, hi %f' % (lo, mx, hi)) 64 warn = 'WARNING: sky estimate not bracketed by peak: lo %f, sky %f, hi %f' % (lo, mx, hi) 65 66 if return_fit: 67 bfit = X[0] + X[1] * x + X[2] * x**2 68 return (x * xscale + x0, b, bfit, mx, warn, binedges1,counts1) 69 70 return mx 71 72 73def parse_ranges(s): 74 ''' 75 Parse PBS job array-style ranges: NNN,MMM-NNN,PPP 76 77 *s*: string 78 79 Returns: [ int, int, ... ] 80 ''' 81 tiles = [] 82 words = s.split() 83 for w in words: 84 for a in w.split(','): 85 if '-' in a: 86 aa = a.split('-') 87 if len(aa) != 2: 88 raise RuntimeError('With an arg containing a dash, expect two parts, in word "%s"' % a) 89 start = int(aa[0]) 90 end = int(aa[1]) 91 for i in range(start, end+1): 92 tiles.append(i) 93 else: 94 tiles.append(int(a)) 95 return tiles 96 97 98def patch_image(img, mask, dxdy = [(-1,0),(1,0),(0,-1),(0,1)], 99 required=None): 100 ''' 101 Patch masked pixels by iteratively averaging non-masked neighboring pixels. 102 103 WARNING: this modifies BOTH the "img" and "mask" arrays! 104 105 mask: True for good pixels 106 required: if non-None: True for pixels you want to be patched. 107 dxdy: Pixels to average in, relative to pixels to be patched. 108 109 Returns True if patching was successful. 110 ''' 111 assert(img.shape == mask.shape) 112 assert(len(img.shape) == 2) 113 h,w = img.shape 114 Nlast = -1 115 while True: 116 needpatching = np.logical_not(mask) 117 if required is not None: 118 needpatching *= required 119 I = np.flatnonzero(needpatching) 120 if len(I) == 0: 121 break 122 if len(I) == Nlast: 123 return False 124 #print 'Patching', len(I), 'pixels' 125 Nlast = len(I) 126 iy,ix = np.unravel_index(I, img.shape) 127 psum = np.zeros(len(I), img.dtype) 128 pn = np.zeros(len(I), int) 129 130 for dx,dy in dxdy: 131 ok = True 132 if dx < 0: 133 ok = ok * (ix >= (-dx)) 134 if dx > 0: 135 ok = ok * (ix <= (w-1-dx)) 136 if dy < 0: 137 ok = ok * (iy >= (-dy)) 138 if dy > 0: 139 ok = ok * (iy <= (h-1-dy)) 140 141 # darn, NaN * False = NaN, not zero. 142 finite = np.isfinite(img [iy[ok]+dy, ix[ok]+dx]) 143 ok[ok] *= finite 144 145 psum[ok] += (img [iy[ok]+dy, ix[ok]+dx] * 146 mask[iy[ok]+dy, ix[ok]+dx]) 147 pn[ok] += mask[iy[ok]+dy, ix[ok]+dx] 148 149 # print 'ix', ix 150 # print 'iy', iy 151 # print 'dx,dy', dx,dy 152 # print 'ok', ok 153 # print 'psum', psum 154 # print 'pn', pn 155 156 img.flat[I] = (psum / np.maximum(pn, 1)).astype(img.dtype) 157 mask.flat[I] = (pn > 0) 158 #print 'Patched', np.sum(pn > 0) 159 return True 160 161def clip_wcs(wcs1, wcs2, makeConvex=True, pix1=None, pix2=None): 162 ''' 163 Returns a pixel-space polygon in WCS1 after it is clipped by WCS2. 164 Returns an empty list if the two WCS headers do not intersect. 165 166 Note that due to weakness in the clip_polygon method, wcs2 must be convex. 167 If makeConvex=True, we find the convex hull and clip with that. 168 169 If pix1 is not None, pix1=(xx,yy), xx and yy lists of pixel coordinates defining 170 the boundary of the image, in CLOCKWISE order; default is the edges and midpoints. 171 172 ''' 173 if pix1 is None: 174 w1,h1 = wcs1.get_width(), wcs1.get_height() 175 x1,x2,x3 = 0.5, w1/2., w1+0.5 176 y1,y2,y3 = 0.5, h1/2., h1+0.5 177 xx = [x1, x1, x1, x2, x3, x3, x3, x2] 178 yy = [y1, y2, y3, y3, y3, y2, y1, y1] 179 else: 180 xx,yy = pix1 181 #rr,dd = wcs1.pixelxy2radec(xx, yy) 182 183 if pix2 is None: 184 w2,h2 = wcs2.get_width(), wcs2.get_height() 185 x1,x2,x3 = 0.5, w2/2., w2+0.5 186 y1,y2,y3 = 0.5, h2/2., h2+0.5 187 XX = [x1, x1, x1, x2, x3, x3, x3, x2] 188 YY = [y1, y2, y3, y3, y3, y2, y1, y1] 189 else: 190 XX,YY = pix2 191 rr,dd = wcs2.pixelxy2radec(XX, YY) 192 ok,XX,YY = wcs1.radec2pixelxy(rr, dd) 193 # XX,YY is the clip polygon in wcs1 pixel coords. 194 # Not necessarily clockwise at this point! 195 196 #print 'XX,YY', XX, YY 197 198 if makeConvex: 199 from scipy.spatial import ConvexHull 200 points = np.vstack((XX,YY)).T 201 #print 'points', points.shape 202 hull = ConvexHull(points) 203 #print 'Convex hull:', hull 204 #print hull.vertices 205 206 # ConvexHull returns the vertices in COUNTER-clockwise order. Reverse. 207 v = np.array(list(reversed(hull.vertices))).astype(int) 208 XX = XX[v] 209 YY = YY[v] 210 211 #plt.plot(XX, YY, 'm-') 212 #plt.plot(XX[0], YY[0], 'mo') 213 #plt.savefig('clip2.png') 214 else: 215 # Ensure points are listed in CLOCKWISE order. 216 crosses = [] 217 for i in range(len(XX)): 218 j = (i + 1) % len(XX) 219 k = (i + 2) % len(XX) 220 xi,yi = XX[i], YY[i] 221 xj,yj = XX[j], YY[j] 222 xk,yk = XX[k], YY[k] 223 dx1, dy1 = xj - xi, yj - yi 224 dx2, dy2 = xk - xj, yk - yj 225 cross = dx1 * dy2 - dy1 * dx2 226 #print 'cross', cross 227 crosses.append(cross) 228 crosses = np.array(crosses) 229 #print 'cross products', crosses 230 if np.all(crosses >= 0): 231 # Reverse 232 #print 'Reversing wcs2 points' 233 XX = np.array(list(reversed(XX))) 234 YY = np.array(list(reversed(YY))) 235 236 clip = clip_polygon(list(zip(xx, yy)), list(zip(XX, YY))) 237 clip = np.array(clip) 238 239 if False: 240 import pylab as plt 241 plt.clf() 242 plt.plot(xx, yy, 'b.-') 243 plt.plot(xx[0], yy[0], 'bo') 244 plt.plot(XX, YY, 'r.-') 245 plt.plot(XX[0], YY[0], 'ro') 246 if len(clip) > 0: 247 plt.plot(clip[:,0], clip[:,1], 'm.-', alpha=0.5, lw=2) 248 plt.savefig('clip1.png') 249 250 return clip 251 252 253 254def polygon_area(poly): 255 ''' 256 NOTE, unlike many of the other methods in this module, takes: 257 258 poly = (xx,yy) 259 260 where xx,yy MUST repeat the starting point at the end of the polygon. 261 ''' 262 xx,yy = poly 263 x,y = np.mean(xx), np.mean(yy) 264 area = 0. 265 for dx0,dy0,dx1,dy1 in zip(xx-x, yy-y, xx[1:]-x, yy[1:]-y): 266 # area: 1/2 cross product 267 area += np.abs(dx0 * dy1 - dx1 * dy0) 268 return 0.5 * area 269 270def clip_polygon(poly1, poly2): 271 ''' 272 Returns a new polygon resulting from taking poly1 and clipping it 273 to lie inside poly2. 274 275 WARNING, the polygons must be listed in CLOCKWISE order. 276 277 WARNING, the clipping polygon, poly2, must be CONVEX. 278 ''' 279 # from clipper import Clipper, Point, PolyType, ClipType, PolyFillType 280 # ''' 281 # ''' 282 # c = Clipper() 283 # p1 = [Point(x,y) for x,y in poly1] 284 # p2 = [Point(x,y) for x,y in poly2] 285 # c.AddPolygon(p1, PolyType.Subject) 286 # c.AddPolygon(p2, PolyType.Clip) 287 # solution = [] 288 # pft = PolyFillType.EvenOdd 289 # result = c.Execute(ClipType.Intersection, solution, pft, pft) 290 # if len(solution) > 1: 291 # raise RuntimeError('Polygon clipping results in non-simple polygon') 292 # assert(result) 293 # #print 'Result:', result 294 # #print 'Solution:', solution 295 # return [(s.x, s.y) for s in solution[0]] 296 297 # Sutherland-Hodgman algorithm -- thanks, Wikipedia! 298 N2 = len(poly2) 299 # clip by each edge in turn. 300 for j in range(N2): 301 # target "left_right" value 302 clip1 = poly2[j] 303 clip2 = poly2[(j+1)%N2] 304 LRinside = _left_right(clip1, clip2, poly2[(j+2)%N2]) 305 # are poly vertices inside or outside the clip polygon? 306 isinside = [_left_right(clip1, clip2, p) == LRinside 307 for p in poly1] 308 # the resulting clipped polygon 309 clipped = [] 310 N1 = len(poly1) 311 for i in range(N1): 312 S = poly1[i] 313 E = poly1[(i+1)%N1] 314 Sin = isinside[i] 315 Ein = isinside[(i+1)%N1] 316 if Ein: 317 if not Sin: 318 clipped.append(line_intersection(clip1, clip2, S, E)) 319 clipped.append(E) 320 else: 321 if Sin: 322 clipped.append(line_intersection(clip1, clip2, S, E)) 323 poly1 = clipped 324 return poly1 325 326 327def polygons_intersect(poly1, poly2): 328 ''' 329 Determines whether the given 2-D polygons intersect. 330 331 poly1, poly2: np arrays with shape (N,2) 332 ''' 333 334 # Check whether any points in poly1 are inside poly2, 335 # or vice versa. 336 for (px,py) in poly1: 337 if point_in_poly(px,py, poly2): 338 return (px,py) 339 for (px,py) in poly2: 340 if point_in_poly(px,py, poly1): 341 return (px,py) 342 343 # Check for intersections between line segments. O(n^2) brutish 344 N1 = len(poly1) 345 N2 = len(poly2) 346 347 for i in range(N1): 348 for j in range(N2): 349 xy = line_segments_intersect(poly1[i % N1, :], poly1[(i+1) % N1, :], 350 poly2[j % N2, :], poly2[(j+1) % N2, :]) 351 if xy: 352 return xy 353 return False 354 355 356def line_segments_intersect(xy1, xy2, xy3, xy4): 357 ''' 358 Determines whether the two given line segments intersect; 359 360 (x1,y1) to (x2,y2) 361 and 362 (x3,y3) to (x4,y4) 363 ''' 364 (x1,y1) = xy1 365 (x2,y2) = xy2 366 (x3,y3) = xy3 367 (x4,y4) = xy4 368 x,y = line_intersection((x1,y1),(x2,y2),(x3,y3),(x4,y4)) 369 if x is None: 370 # Parallel lines 371 return False 372 if x1 == x2: 373 p1,p2 = y1,y2 374 p = y 375 else: 376 p1,p2 = x1,x2 377 p = x 378 379 if not ((p >= min(p1,p2)) and (p <= max(p1,p2))): 380 return False 381 382 if x3 == x4: 383 p1,p2 = y3,y4 384 p = y 385 else: 386 p1,p2 = x3,x4 387 p = x 388 389 if not ((p >= min(p1,p2)) and (p <= max(p1,p2))): 390 return False 391 return (x,y) 392 393 394def line_intersection(xy1, xy2, xy3, xy4): 395 ''' 396 Determines the point where the lines described by 397 (x1,y1) to (x2,y2) 398 and 399 (x3,y3) to (x4,y4) 400 intersect. 401 402 Note that this may be beyond the endpoints of the line segments. 403 404 Probably raises an exception if the lines are parallel, or does 405 something numerically crazy. 406 ''' 407 (x1,y1) = xy1 408 (x2,y2) = xy2 409 (x3,y3) = xy3 410 (x4,y4) = xy4 411 # This code started with the equation from Wikipedia, 412 # then I added special-case handling. 413 # bottom = ((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)) 414 # if bottom == 0: 415 # raise RuntimeError("divide by zero") 416 # t1 = (x1 * y2 - y1 * x2) 417 # t2 = (x3 * y4 - y3 * x4) 418 # px = (t1 * (x3 - x4) - t2 * (x1 - x2)) / bottom 419 # py = (t1 * (y3 - y4) - t2 * (y1 - y2)) / bottom 420 421 # From http://wiki.processing.org/w/Line-Line_intersection 422 bx = float(x2) - float(x1) 423 by = float(y2) - float(y1) 424 dx = float(x4) - float(x3) 425 dy = float(y4) - float(y3) 426 b_dot_d_perp = bx*dy - by*dx 427 if b_dot_d_perp == 0: 428 return None,None 429 cx = float(x3) - float(x1) 430 cy = float(y3) - float(y1) 431 t = (cx*dy - cy*dx) / b_dot_d_perp 432 return x1 + t*bx, y1 + t*by 433 434def _left_right(xy1, xy2, xy3): 435 ''' 436 is (x3,y3) to the 'left' or 'right' of the line from (x1,y1) to (x2,y2) ? 437 ''' 438 (x1,y1) = xy1 439 (x2,y2) = xy2 440 (x3,y3) = xy3 441 dx2,dy2 = x2-x1, y2-y1 442 dx3,dy3 = x3-x1, y3-y1 443 return (dx2 * dy3 - dx3 * dy2) > 0 444 445 446def point_in_poly(x, y, poly): 447 ''' 448 Performs a point-in-polygon test for numpy arrays of *x* and *y* 449 values, and a polygon described as 2-d numpy array (with shape (N,2)) 450 451 poly: N x 2 array 452 453 Returns a numpy array of bools. 454 ''' 455 x = np.atleast_1d(x) 456 y = np.atleast_1d(y) 457 inside = np.zeros(x.shape, bool) 458 # This does a winding test -- count how many times a horizontal ray 459 # from (-inf,y) to (x,y) crosses the boundary. 460 for i in range(len(poly)): 461 j = (i-1 + len(poly)) % len(poly) 462 xi,xj = poly[i,0], poly[j,0] 463 yi,yj = poly[i,1], poly[j,1] 464 465 if yi == yj: 466 continue 467 468 I = np.logical_and( 469 np.logical_or(np.logical_and(yi <= y, y < yj), 470 np.logical_and(yj <= y, y < yi)), 471 x < (xi + ((xj - xi) * (y - yi) / (yj - yi)))) 472 inside[I] = np.logical_not(inside[I]) 473 return inside 474 475def lanczos_filter(order, x, out=None): 476 x = np.atleast_1d(x) 477 nz = np.logical_and(x != 0., np.logical_and(x < order, x > -order)) 478 nz = np.flatnonzero(nz) 479 if out is None: 480 out = np.zeros(x.shape, dtype=float) 481 else: 482 out[x <= -order] = 0. 483 out[x >= order] = 0. 484 pinz = pi * x.flat[nz] 485 out.flat[nz] = order * np.sin(pinz) * np.sin(pinz / order) / (pinz**2) 486 out[x == 0] = 1. 487 return out 488 489def get_overlapping_region(xlo, xhi, xmin, xmax): 490 ''' 491 Given a range of integer coordinates that you want to, eg, cut out 492 of an image, [xlo, xhi], and bounds for the image [xmin, xmax], 493 returns the range of coordinates that are in-bounds, and the 494 corresponding region within the desired cutout. 495 496 For example, say you have an image of shape H,W and you want to 497 cut out a region of halfsize "hs" around pixel coordinate x,y, but 498 so that coordinate x,y is centered in the cutout even if x,y is 499 close to the edge. You can do: 500 501 cutout = np.zeros((hs*2+1, hs*2+1), img.dtype) 502 iny,outy = get_overlapping_region(y-hs, y+hs, 0, H-1) 503 inx,outx = get_overlapping_region(x-hs, x+hs, 0, W-1) 504 cutout[outy,outx] = img[iny,inx] 505 506 ''' 507 if xlo > xmax or xhi < xmin or xlo > xhi or xmin > xmax: 508 return ([], []) 509 510 assert(xlo <= xhi) 511 assert(xmin <= xmax) 512 513 xloclamp = max(xlo, xmin) 514 Xlo = xloclamp - xlo 515 516 xhiclamp = min(xhi, xmax) 517 Xhi = Xlo + (xhiclamp - xloclamp) 518 519 #print 'xlo, xloclamp, xhiclamp, xhi', xlo, xloclamp, xhiclamp, xhi 520 assert(xloclamp >= xlo) 521 assert(xloclamp >= xmin) 522 assert(xloclamp <= xmax) 523 assert(xhiclamp <= xhi) 524 assert(xhiclamp >= xmin) 525 assert(xhiclamp <= xmax) 526 #print 'Xlo, Xhi, (xmax-xmin)', Xlo, Xhi, xmax-xmin 527 assert(Xlo >= 0) 528 assert(Xhi >= 0) 529 assert(Xlo <= (xhi-xlo)) 530 assert(Xhi <= (xhi-xlo)) 531 532 return (slice(xloclamp, xhiclamp+1), slice(Xlo, Xhi+1)) 533 534 535 536if __name__ == '__main__': 537 import matplotlib 538 matplotlib.use('Agg') 539 import pylab as plt 540 import numpy as np 541 from astrometry.util.plotutils import * 542 ps = PlotSequence('miscutils') 543 544 np.random.seed(42) 545 546 if True: 547 p2 = np.array([[0,0],[0,4],[4,4],[4,0]]) 548 549 for i,p1 in enumerate([ np.array([[0,0],[0,2],[2,2],[2,0]]), 550 np.array([[-1,-1],[0,2],[2,2],[2,0]]), 551 np.array([[4,0],[0,4],[-4,0],[0,-4]]), 552 np.array([[-1,2],[2,5],[5,2],[2,-1]]), 553 ] + [None]*10): 554 if p1 is None: 555 p1 = np.random.uniform(high=6., low=-2, size=(4,2)) 556 pc = np.array(clip_polygon(p1, p2)) 557 plt.clf() 558 I = np.array([0,1,2,3,0]) 559 plt.plot(p1[I,0], p1[I,1], 'b-', lw=3, alpha=0.5) 560 plt.plot(p2[I,0], p2[I,1], 'k-') 561 I = np.array(range(len(pc)) + [0]) 562 plt.plot(pc[I,0], pc[I,1], 'r-') 563 plt.axis([-1,5,-1,5]) 564 plt.savefig('clip-%02i.png' % i) 565 import sys 566 sys.exit(0) 567 568 if True: 569 for i in range(20): 570 if i <= 10: 571 xy1 = np.array([[0,0],[0,4],[4,4],[4,0]]) 572 else: 573 xy1 = np.random.uniform(high=10., size=(4,2)) 574 xy2 = np.random.uniform(high=10., size=(4,2)) 575 plt.clf() 576 I = np.array([0,1,2,3,0]) 577 xy = polygons_intersect(xy1, xy2) 578 if xy: 579 cc = 'r' 580 x,y = xy 581 plt.plot(x,y, 'mo', mec='m', mfc='none', ms=20, mew=3, zorder=30) 582 else: 583 cc = 'k' 584 plt.plot(xy1[I,0], xy1[I,1], '-', color=cc, zorder=20, lw=3) 585 plt.plot(xy2[I,0], xy2[I,1], '-', color=cc, zorder=20, lw=3) 586 ax = plt.axis() 587 plt.axis([ax[0]-0.5, ax[1]+0.5, ax[2]-0.5, ax[3]+0.5]) 588 ps.savefig() 589 590 if False: 591 X,Y = np.meshgrid(np.linspace(-1,11, 20), np.linspace(-1,11, 23)) 592 X = X.ravel() 593 Y = Y.ravel() 594 for i in range(20): 595 if i == 0: 596 xy = np.array([[0,0],[0,10],[10,10],[10,0]]) 597 else: 598 xy = np.random.uniform(high=10., size=(4,2)) 599 plt.clf() 600 I = np.array([0,1,2,3,0]) 601 plt.plot(xy[I,0], xy[I,1], 'r-', zorder=20, lw=3) 602 inside = point_in_poly(X, Y, xy) 603 plt.plot(X[inside], Y[inside], 'bo') 604 out = np.logical_not(inside) 605 plt.plot(X[out], Y[out], 'ro') 606 ax = plt.axis() 607 plt.axis([ax[0]-0.5, ax[1]+0.5, ax[2]-0.5, ax[3]+0.5]) 608 ps.savefig() 609 610 611 612 if True: 613 # intersection() 614 for i in range(20): 615 if i == 0: 616 x1 = x2 = 0 617 y1 = 0 618 y2 = 1 619 x3 = 1 620 x4 = -1 621 y3 = 0 622 y4 = 1 623 elif i == 1: 624 x1,y1 = 0,0 625 x2,y2 = 0,1 626 x3,y3 = -3,0 627 x4,y4 = -2,0 628 elif i == 2: 629 x1,y1 = 1,0 630 x2,y2 = 0,1 631 x3,y3 = -3,0 632 x4,y4 = -2,0 633 elif i == 3: 634 x1,y1 = 0,1 635 x2,y2 = 1,0 636 x3,y3 = 0,-3 637 x4,y4 = 0,-2 638 elif i == 4: 639 x1,y1 = 0,0 640 x2,y2 = 0,1 641 x3,y3 = 0,2 642 x4,y4 = 0,3 643 elif i == 5: 644 x1,y1 = -1,0 645 x2,y2 = 1, 0 646 x3,y3 = 0, 2 647 x4,y4 = 0.5, 1 648 else: 649 xy = np.random.uniform(high=10., size=(8,)) 650 x1,y1,x2,y2,x3,y3,x4,y4 = xy 651 plt.clf() 652 plt.plot([x1,x2],[y1,y2], 'r-', zorder=20, lw=3) 653 plt.plot([x3,x4],[y3,y4], 'b-', zorder=20, lw=3) 654 x,y = line_intersection((x1,y1),(x2,y2),(x3,y3),(x4,y4)) 655 plt.plot(x, y, 'kx', ms=20, zorder=25) 656 plt.plot([x1,x],[y1,y], 'k--', alpha=0.5, zorder=15) 657 plt.plot([x2,x],[y2,y], 'k--', alpha=0.5, zorder=15) 658 plt.plot([x3,x],[y3,y], 'k--', alpha=0.5, zorder=15) 659 plt.plot([x4,x],[y4,y], 'k--', alpha=0.5, zorder=15) 660 661 # line_segments_intersect() 662 if line_segments_intersect((x1,y1),(x2,y2),(x3,y3),(x4,y4)): 663 plt.plot(x,y, 'mo', mec='m', mfc='none', ms=20, mew=3, zorder=30) 664 ax = plt.axis() 665 plt.axis([ax[0]-0.5, ax[1]+0.5, ax[2]-0.5, ax[3]+0.5]) 666 ps.savefig() 667