""" Module to create tessellations This algorithm successfully crops the diagram and is capable of dealing with vertical lines (8 June 2018) """ import warnings import numpy as np import scipy.spatial as sp import triangle as tgl from collections import OrderedDict import matplotlib.patches as mp from matplotlib.collections import PatchCollection from matplotlib.ticker import LogFormatter from matplotlib.colors import LogNorm, Normalize import matplotlib.pyplot as plt import algebra as alg import geometry as geo import tools import workspace as ws # Ignore warnings warnings.filterwarnings('ignore') class PowerDiagram(): """ Class for a Power diagram, or a Laguerre-Voronoi diagram """ def __init__(self, x, y, r, contour=None): self.x = x self.y = y self.r = r self.circ_areas = np.pi * r ** 2 self.nc = x.size self.contour = contour self.any_errors = False # self.build() # self.get_polygons() def build(self): """ Build it """ ws.log("Building the power diagram") x, y, r = self.x, self.y, self.r # Generating circles gc = GCircles(x, y, r) gc.triangulate() # First iteration pairs, trios, segments, pairs_trios, pairs_oprs, pairs_opts, \ pairs_allpts, pvpts, lines = build_laguerre(gc, gc.trios) # First correction new_trios, new_pairs, rmv_trios, rmv_pairs, finished = \ one_correction(pairs_oprs, pairs_opts, pairs_trios, segments) trios = update_trios(new_trios, rmv_trios, trios) # Iterate until there's no more bad crossings # Iteration counter corrcount = 0 while not finished: pairs, trios, segments, pairs_trios, pairs_oprs, pairs_opts, \ pairs_allpts, pvpts, lines = build_laguerre(gc, trios) new_trios, new_pairs, rmv_trios, rmv_pairs, finished = \ one_correction(pairs_oprs, pairs_opts, pairs_trios, segments) trios = update_trios(new_trios, rmv_trios, trios) # Go fix the only-two-connections issue poss_lost = [] poss_lost, rmv_trios, new_trios = fix_biconnections(poss_lost, pairs, gc, pairs_trios) if len(poss_lost) > 0: trios = update_trios(new_trios, rmv_trios, trios) pairs, trios, segments, pairs_trios, pairs_oprs, pairs_opts, \ pairs_allpts, pvpts, lines = build_laguerre(gc, trios) # 'This should be fixed now' whoslost = poss_lost[0] corrcount += 1 enough = corrcount == gc.nc / 4 # Check if this is finished finished = finished or enough ws.log('Number of corrections: %i'%corrcount) if enough: error_msg = 'TESSELLATION ERROR: Something went ' \ + 'wrong with this set of circles\n'\ + 'This is probably due to a bug in the algorithm\n' \ + 'Please try a different set of circles\n' ws.log(error_msg) self.any_errors = True print('TESSELLATION ERROR') return ws.log('Correction iterations: %i'%corrcount) # Check if two connections cross for ip, pair1 in enumerate(pairs): a = segments[pair1] for pair2 in pairs[ip+1:]: b = segments[pair2] if isinstance(a, geo.Segment) and isinstance(b, geo.Segment): crossing, _ = geo.crossing_segments(a, b) overlapping = geo.overlapping_segments(a, b) if crossing or overlapping: ws.log('TESSELLATION ERROR: An error was found for pairs: %s, %s'%(str(pair1), str(pair2))) ws.log('This is due to a bug in the algorithm') ws.log('Please try a different set of circles') self.any_errors = True print('TESSELLATION ERROR') return else: pass # If finished, crop the diagram using the contour if self.contour is not None: # print('') # print('segments:') # for k, v in segments.items(): # print(k, v) # Crop it segments = crop_diagram(gc, pairs, segments, pvpts, lines, self.contour) # Finally, remove the segments which are not geo.Segment # instances for k, seg in segments.items(): if not isinstance(seg, geo.Segment): if len(seg) != 2: segments[k] = None if len(seg) == 2: segments[k] = geo.Segment(seg) self.pairs = pairs self.segments = segments self.trios = trios self.get_polygons() def build_preexisting(self, pairs, segments): """ Set up a diagram with preset properties Created on 4 October 2018 """ pairs = list(pairs) self.pairs = pairs self.segments = segments trios = [] for i, p1 in enumerate(pairs): for p2 in pairs[i+1:]: s1 = set(p1) s2 = set(p2) # If there is only three points involved between the two # pairs, this is clearly a potential trio union = s1.union(s2) if len(union) == 3: # But not all trios are real in our diagram, so we need to # check that this trio exists really # For this, check that the two non-common points, e.g. # (2, 3) if the pairs were (1, 2) and (1, 3) are a listed # pair if tuple(s1.symmetric_difference(s2)) in pairs: trios.append(tuple(union)) self.trios = trios self.get_polygons() def get_polygons(self): """ Create the polygons for each cell """ segments = self.segments x, y = self.x, self.y nc = self.nc cellverts = OrderedDict() polygons = OrderedDict() polg_areas = np.zeros(nc) # Add vertices to the cells for (i, j), seg in segments.items(): try: a, b = seg.a, seg.b except AttributeError: # This means that some segment was [None] print("None segment for (%i, %i)"%(i, j)) pass else: cellverts = tools.append_items(cellverts, i, [a, b]) cellverts = tools.append_items(cellverts, j, [a, b]) # Set vertices for i in range(nc): points = np.array(cellverts[i]) # Reshape if points.shape[-1] == 1: points = points.reshape(points.shape[:-1]) # Set them points = np.unique(points, axis=0) # Sort points by angle cellpoint = (x[i], y[i]) # Sort them by angle. Use some averaged centroid rather than # the cellpoint itself. Because if cellpoint is included, it # may be confusing and the sorting can go wrong # This is OK, since the polygon is always convex points = geo.sort_by_angle(points) cellverts[i] = points # Create polygon pol = geo.Polygon(points) # Check that the polygon contains the point or this is on # the contour, at least. If not, add it to the polygon polccell = pol.contains_point(cellpoint) # if False: if not polccell: points = np.append(points.tolist(), [list(cellpoint)], axis=0) points = geo.sort_by_angle(points) cellverts[i] = points pol = geo.Polygon(points) try: polg_areas[i] = pol.area except AttributeError: # This is not a polygon pass else: # If it is, save it polygons[i] = pol # Store in attributes self.cellverts = cellverts self.polygons = polygons self.polg_areas = polg_areas self.free_areas = polg_areas - self.circ_areas def draw( self, ax, draw_polygons=True, colour='k', line_alpha=1., fill_alpha=1., linestyle='-', linewidth=1., markers=None, values=None, precompv=None, vmin=None, vmax=None, title=None, cmap=None, cblabel=None, logscale=False, allowed_cbticklabels=(None,), extend='neither', zorder=0 ): """ Draw the diagram """ drawvalues = values is not None if drawvalues: patches = [] colours = np.zeros(self.nc) for i, pol in self.polygons.items(): # Colours and filled polygons if drawvalues: filled_pol = mp.Polygon(self.cellverts[i]) patches.append(filled_pol) if values == 'areas': colours[i] = pol.area if values == 'inv_areas': colours[i] = 1. / self.free_areas[i] if values == 'free_areas': colours[i] = self.free_areas[i] if values == 'precomputed': if logscale: colours[i] = np.abs(precompv[i]) else: colours[i] = precompv[i] # Colour for the polygon contours if colour == 'same_as_facecolors': pc_ = colours[i] else: pc_ = colour # Draw the polygon if draw_polygons: pol.draw(ax=ax, colour=pc_, alpha=line_alpha, linestyle=linestyle, linewidth=linewidth, markers=markers, zorder=zorder) # Line colours if drawvalues: if vmin is None: vmin = colours.min() if vmax is None: vmax = colours.max() if logscale: vmin = np.abs(vmin) vmax = np.abs(vmax) norm = LogNorm(vmin=vmin, vmax=vmax) else: norm = Normalize(vmin=vmin, vmax=vmax) m = plt.cm.ScalarMappable(norm=norm, cmap=cmap) line_colours = m.to_rgba(colours) print(line_colours) if drawvalues: collection = PatchCollection(patches, cmap=cmap, norm=norm) ax.add_collection(collection) collection.set_alpha(fill_alpha) collection.set_array(colours) if logscale: collection.set_clim(vmin=vmin, vmax=vmax) else: collection.set_clim(vmin=0., vmax=vmax) if colour == 'same_as_facecolors': collection.set_edgecolors(line_colours) else: collection.set_edgecolor(colour) collection.set_linewidth(linewidth) if values == 'precomputed': collection.set_clim(vmin=vmin, vmax=vmax) # Color bar if logscale: l_f = LogFormatter(10, labelOnlyBase=False) cb = plt.colorbar(collection, ax=ax, orientation='vertical', format=l_f, extend=extend) # Set minor ticks # We need to nomalize the tick locations so that they're in the range from 0-1... # minorticks = p.norm(np.arange(1, 10, 2)) # Minimum and maximum exponents min_exp = int(np.log10(min(vmin, vmax))) max_exp = int(np.log10(max(vmin, vmax))) + 1 # Ticks for the colorbar in case it's logscale minorticks = 10. ** np.arange(min_exp, max_exp + 1, 1) minorticks_ = minorticks.tolist() for i in range(2, 10, 1): minorticks_ = minorticks_ + (float(i) * minorticks).tolist() minorticks_.sort() minorticks = np.array(minorticks_) minorticks = minorticks[np.where(minorticks >= vmin)] minorticks = minorticks[np.where(minorticks <= vmax)] print(minorticks) cb.set_ticks(minorticks) cb.set_ticklabels([str(int(x)) if x in allowed_cbticklabels else '' for x in minorticks]) else: cb = plt.colorbar(collection, ax=ax, orientation='vertical', extend=extend) cb.set_label(cblabel) if title is not None: ax.set_title(title) ax.set_xlabel(r'$\rm x$ ($\rm \mu m$)') ax.set_ylabel(r'$\rm y$ ($\rm \mu m$)') class PointsEnsemble(): """ Cloud of points """ def __init__(self, x, y): self.x = x self.y = y self.points = np.array([self.x, self.y]).T self.nc = len(self.x) # Matrices self.dx = alg.compl_mat(x) self.dy = alg.compl_mat(y) def triangulate(self): """ Get the Delaunay triangulation of the points """ # Triangulation if self.nc == 2: self.pairs = np.array([np.arange(2)]) self.trios = np.array([np.arange(2)]) self.triangles = {} else: tri = sp.Delaunay(self.points) trios = tri.simplices.copy() self.triangulation = tri self.update_triangulation(trios) def update_triangulation(self, trios): """ Update the triangulation for a new arrangement of trios """ pairs = [] # Triangles triangles = {} for trio in trios: trio = tuple(sorted(trio)) triangle = geo.Triangle(self.points, indices=trio, fas=None) triangles[trio] = triangle # Compute the resistors for each side of the triangle for i, pair in enumerate(triangle.pairs): # We do not want to repeat pairs, right? pair = tuple(sorted(pair)) if pair not in pairs: pairs.append(pair) pairs = np.array(pairs) self.pairs = pairs self.triangles = triangles trios = np.sort(trios) self.trios = trios class GCircles(PointsEnsemble): """ Ensemble of generating circles in the system """ def __init__(self, x, y, r): PointsEnsemble.__init__(self, x, y) self.r = r self.nc = self.x.size # Vectors # Sigmas self.sigma = r * r - (x * x + y * y) # Matrices self.ds = alg.compl_mat(self.sigma) class LinesEnsemble(): """ Ensemble of lines """ def __init__(self): pass def cutoffpairs(self): """ Reduce the size of the matrices according to the pairing in the triangulation """ pi, pj = self.pairs.T self.m = self.m[pi, pj] self.n = self.n[pi, pj] self.pi, self.pj = pi, pj def vertical_to_mbform(self): """ Get all the lines in point-slope form, so I can deal with vertical lines in other calculations """ # Provisional b-points (for the point-slope form) self.bx = np.zeros_like(self.n) self.by = self.n # Do the necessary with vertical lines pi, pj = self.pi, self.pj dxvec = self.pe.dx[pi, pj] dyvec = self.pe.dy[pi, pj] dsvec = self.pe.ds[pi, pj] ivlines = np.where(dyvec == 0.) # Find the point bx bx = -0.5 * dsvec[ivlines] / dxvec[ivlines] # by by = self.pe.y[self.pairs[ivlines, 0]] self.bx[ivlines] = bx self.by[ivlines] = by def dict(self): """ Create line instances in the geometry module so that they can be drawn on a figure """ bx, by = self.bx, self.by pairs, n, m = self.pairs, self.n, self.m lines = {} for i, pair in enumerate(pairs): bxi, byi = bx[i], by[i] lines[tuple(pair)] = geo.StraightLine((bxi, byi), m[i]) # Store them as attributes self.lines = lines self.bxvec = bx self.byvec = by self.bpoints = np.array([bx, by]).T class LinesConnectingCenters(LinesEnsemble): """ Ensemble of lines connecting the centers of an ensemble of circles """ def __init__(self, gc): self.pe = gc self.nc = self.pe.nc self.pairs = self.pe.pairs self.nl = len(self.pairs) LinesEnsemble.__init__(self) def build(self): """ Build the lines """ x, y = self.pe.x, self.pe.y dx, dy = self.pe.dx, self.pe.dy m = dy / dx n = alg.vectomat(y).T - m * alg.vectomat(x).T # Store everything as attributes self.m = m self.n = n self.cutoffpairs() # Deal with vertical lines self.vertical_to_mbform() # Points on the lines # Indices of (the first of) the points each line belongs to which = self.pairs[:, 0] self.bx = x[which] self.by = y[which] # Create its instances in the geometry module self.dict() class PowerLines(LinesEnsemble): """ Ensemble of power lines between circles """ def __init__(self, gc, lcc): self.pe = gc self.lcc = lcc self.nc = gc.nc self.pairs = self.pe.pairs self.nl = len(self.pairs) LinesEnsemble.__init__(self) def build(self): """ Build the lines """ dx, dy, ds = self.pe.dx, self.pe.dy, self.pe.ds # Slopes m = - dx / dy # Intercept (y value) n = -0.5 * ds / dy # Store everything as attributes self.m = m self.n = n self.cutoffpairs() # Deal with vertical lines self.vertical_to_mbform() # Create its instances in the geometry module self.dict() # Find points in the lines for geo.StraightLine # Find the power-line mid-points between pairs mps = geo.intersections_many_lines([self, self.lcc]) xm, ym = mps # Only interested in the points that involve both ensembles xm = xm[:self.nl, self.nl:] ym = ym[:self.nl, self.nl:] xm = np.diag(xm) ym = np.diag(ym) self.bx = xm self.by = ym # Create its instances in the geometry module self.dict() def remove_repeated_points(x, y): """ Make a set of the points so there are no repeated points """ # Use my geometry module, because it may very well be that points # which are the same in theory, are not the same in practice # Arrange points points = np.array([x, y]).T # Distances among all points # dists = np.hypot(alg.compl_mat(x), alg.compl_mat(y)) # pairs, dists = alg.one_side(dists) pairs, dists = geo.dist_cloud(x, y) # Where are the points nearly identical where = np.where(dists <= 1.e-9) pairs_matched = pairs[where] # Remove the second one from each pair of nearly identical points allpts = set(pairs.flatten().tolist()) allmatched = set(pairs_matched.flatten().tolist()) removed = set(pairs_matched[:, 1].tolist()) survivors = allpts - removed # survivors = tuple(survivors) survivors = np.array(list(survivors)) points = points[survivors] x, y = points.T # If there is only one point, make x and y be arrays if points.ndim == 1: x = np.array([x]) y = np.array([y]) return x, y def build_laguerre(gc, trios): """ Build the Laguerre diagram from a cloud of points and its given triangulation """ # Update the triangulation gc.update_triangulation(trios) trios = gc.trios # Lines connecting centers lcc = LinesConnectingCenters(gc) lcc.build() # Power lines pl = PowerLines(gc, lcc) pl.build() # Possible Vertices (PV) # Intersections among power lines xi, yi = geo.intersections_many_lines([pl]) lpairs, xi = alg.one_side(xi) lpairs, yi = alg.one_side(yi) # Determine the group of circles for each possible vertex pvnc = len(yi) gcpvc = pl.pairs[lpairs] gcpvc = gcpvc.reshape((pvnc, 4)) gcpvc = np.sort(gcpvc) lgp = np.array(list(map(lambda l: len(set(l)), gcpvc))) keep_these = np.where(lgp == 3) xi = xi[keep_these] yi = yi[keep_these] lpairs = lpairs[keep_these] gcpvc = gcpvc[keep_these] pvnc = len(xi) gcpvc = np.array(list(map(lambda l: sorted(list(set(l))), gcpvc))) # Now I have to iterate over trios # Coordinates of the possible vertex for each triangle xpv = {} ypv = {} # Segment for each pair segments = {} # Number of PV points for each pair pvpts = {} # One (if the pair is on the boundary) or two trios per pair pairs_trios = {} # Other pairs of the quad whose diagonal is a certain pair pairs_oprs = {} # Other points of the quad whose diagonal is a certain pair pairs_opts = {} # All points of the quad whose diagonal is a certain pair pairs_allpts = {} # Iterate for trio in trios: trio_extended = np.tile(trio, pvnc).reshape(gcpvc.shape) matches = np.equal(gcpvc, trio_extended) matches = np.all(matches, axis=1) matches = np.where(matches == True)[0] matches = np.unique(matches) xim = xi[matches] yim = yi[matches] # Now reduce the duplicated points for this triangle if len(xim) > 1: xim, yim = remove_repeated_points(xim, yim) # Store them if len(xim) > 0: tt = tuple(trio) xim, yim = xim.mean(), yim.mean() pairs_sorted = [tuple(sorted(item)) for item in gc.triangles[tt].pairs] # Pairs for pair in pairs_sorted: other_pairs = [item for item in pairs_sorted if item != pair] try: segments[pair].append((xim, yim)) pvpts[pair] += 1 except KeyError: segments[pair] = [(xim, yim)] pvpts[pair] = 1 # For this pair, store the trio and the other pairs pairs_trios[pair] = [tt] pairs_oprs[pair] = other_pairs else: pairs_trios[pair].append(tt) pairs_oprs[pair] = pairs_oprs[pair] + other_pairs pa = set(np.array(pairs_oprs[pair]).flatten().tolist() \ + list(pair)) pairs_allpts[pair] = list(pa) pairs_opts[pair] = sorted(list(pa - set(pair))) # List pairs pairs = list(pairs_oprs.keys()) # for k, v in pvpts.items(): # print(k, v) # Get far points segments = get_far_points(gc, segments, pvpts, pl.lines) # Update segments dictionary for k, seg in segments.items(): # if len(seg) != 2: if not isinstance(seg, geo.Segment): if len(seg) < 2: segments[k] = seg # pass if len(seg) == 2: segments[k] = geo.Segment(seg) # Return return pairs, trios, segments, pairs_trios, pairs_oprs, pairs_opts, pairs_allpts, pvpts, pl.lines def get_far_points(gc, segments, pvpts, lines): """ Get the points being far away """ # Pairs that only have one point (on the boundary) iwhobdr = np.where(np.array(list(pvpts.values())) == 1) whobdr = np.array(list(pvpts.keys()))[iwhobdr] # Give them a second point very far away # Calculate limits xmin, xmax = gc.x.min(), gc.x.max() ymin, ymax = gc.y.min(), gc.y.max() dx = xmax - xmin dy = ymax - ymin xfar_left = xmin - 1e3 * dx xfar_rght = xmax + 1e3 * dx yfar_bot = ymin - 1e3 * dy yfar_top = ymax + 1e3 * dy for pair in whobdr: # Format the pair to amke it a tuple pair = tuple(pair.tolist()) # Get point and clean it try: point = tuple([float(x) for x in segments[pair][0]]) except TypeError: point = segments[pair].a if point is not None: powl = lines[pair] # If it's a vertical line, do it differently slope = powl.slope vline = (slope == np.inf) or (slope == -np.inf) or (np.isnan(slope)) if vline: xfar = point[0] ypnt = point[1] farpt_up = geo.Segment([(xfar, ypnt), (xfar, yfar_top)]) farpt_dn = geo.Segment([(xfar, ypnt), (xfar, yfar_bot)]) possible_segs = [farpt_dn, farpt_up] else: yfar_left = powl.equation(xfar_left) yfar_rght = powl.equation(xfar_rght) # Only one of the segments does not intersect with any other sleft = geo.Segment([point, (xfar_left, yfar_left)]) srght = geo.Segment([point, (xfar_rght, yfar_rght)]) possible_segs = [sleft, srght] # Try the left-hand segment and see if that's the good one valid_segment = np.array([False, False]) for i, tryseg in enumerate(possible_segs): this_one_not_valid = False for othpair, othseg in segments.items(): if (pair != othpair) and isinstance(othseg, geo.Segment): crossing, crosspoint = geo.crossing_segments(tryseg, othseg) overlapping = geo.overlapping_segments(tryseg, othseg) this_one_not_valid = crossing or overlapping if this_one_not_valid: # Not this one break if not this_one_not_valid: valid_segment[i] = True if (valid_segment == True).all(): # Choose the shortest one s0 = possible_segs[0] s1 = possible_segs[1] if s0.length < s1.length: segments[pair] = s0 else: segments[pair] = s1 elif (valid_segment == True).any(): segments[pair] = possible_segs[np.where(valid_segment == True)[0][0]] else: print('No far point could be given to:', pair) return segments def crop_diagram(gc, pairs, segments, pvpts, lines, contour): """ Crop the diagram using a contour """ segments = get_far_points(gc, segments, pvpts, lines) # Calculate limits cont_xy = contour['vertices'] # cont_sg = cont_xy[contour['segments']] contour_polygon = geo.Polygon(cont_xy) x, y = np.array(cont_xy).T xmin, xmax = x.min(), x.max() ymin, ymax = y.min(), y.max() dx = xmax - xmin dy = ymax - ymin # Now get the intersection between a pair's segment and the contour for pair in pairs: seg = segments[pair] if isinstance(seg, geo.Segment): # First off, if the segment is completely outside the contour, # remove it a, b = seg.a, seg.b discard = (not contour_polygon.contains_point(a)) and (not contour_polygon.contains_point(b)) if not discard: for cs in contour_polygon.segments.values(): # xyint = geo.intersection_segments(seg, geo.Segment(cs)) xyint = geo.intersection_segments(seg, cs) if xyint != (np.nan, np.nan): # We found an intersection between the segment and the # contour # Now the point to be replaced is the one falling outside if contour_polygon.contains_point(a): change = 1 elif contour_polygon.contains_point(b): change = 0 # the contour seg.change_point(change, xyint) # We're done with this loop break # Sanity Check after cutting it # If any segment is even partially outside the contour, remove it # I've seen this happen when a segment has one point # outside and one right on the contour, for which no # intersection is found and hence it's left as is # Check if any point is detected to be outside the contour (being # on it may yield True, unfortunately), so in that case, check = (not contour_polygon.contains_point(a)) or (not contour_polygon.contains_point(b)) if check: # roughly check if it's really outside. In that # case, it may be really really far (thus be a # long segment) discard = seg.length > 10. * np.hypot(dx, dy) if discard: ws.log('Discarding %s'%str(pair)) ws.log(str(contour_polygon.contains_point(seg.a))) ws.log(str(seg.points.flatten())) ws.log(str(seg)) ws.log('length = %s, dx = %s, dy = %s'%(str(seg.length), str(dx), str(dy))) ws.log("") segments[pair] = [None] else: # Add the segment segments[pair] = seg else: # Add the segment segments[pair] = seg else: # This is not good doing anymore # segments[pair] = [None] pass # List pairs again pairs = list(segments.keys()) return segments def one_correction(pairs_oprs, pairs_opts, pairs_trios, segments): """ Perform one single trio change """ # New trios new_trios = [] # New pairs new_pairs = [] # Removed trios rmv_trios = [] # Removed pairs rmv_pairs = [] # Iterate over pairs to check crossing segments pairs = list(pairs_trios.keys())[:] for pair in pairs: other_pairs = pairs_oprs[pair] for i, op1 in enumerate(other_pairs): a = segments[op1] for op2 in other_pairs[i+1:]: b = segments[op2] # if (len(a) == 2) & (len(b) == 2): # if True: if isinstance(a, geo.Segment) and isinstance(b, geo.Segment): crossing, _ = geo.crossing_segments(a, b) if crossing: # Change the connection: use the other diagonal of the # quad # Remove bad pair rmv_pairs.append(pair) # Update with new pair new_pair = tuple(pairs_opts[pair]) new_pairs.append(new_pair) # Update its segment xyint = geo.intersection_segments(a, b) try: isgeoSeg = isinstance(segments[new_pair], geo.Segment) except KeyError: segments = tools.append_items(segments, new_pair, [xyint]) else: if not isgeoSeg: segments = tools.append_items(segments, new_pair, [xyint]) # if not isinstance(segments[new_pair], geo.Segment): # try: # segments[new_pair].append(xyint) # except KeyError: # segments[new_pair] = [xyint] # Update the trios # Remove bad trios for trio in pairs_trios[pair]: rmv_trios.append(trio) # Update with new trios for p in pair: new_trio = tuple(sorted(list(new_pair) + [p])) new_trios.append(new_trio) # Update segments dictionary for k, seg in segments.items(): if not isinstance(seg, geo.Segment): # if len(seg) != 2: if len(seg) < 2: segments[k] = seg if len(seg) == 2: segments[k] = geo.Segment(seg) # The important thing here return new_trios, new_pairs, rmv_trios, rmv_pairs, False return new_trios, new_pairs, rmv_trios, rmv_pairs, True def update_trios(new_trios, rmv_trios, trios): """ Update the information according to the correction """ # Set new_trios new_trios = list(set(new_trios)) # Update the lists of pairs and trios trios = [item for item in trios if tuple(item) not in rmv_trios] trios = trios + new_trios # To tuple trios = tuple(map(tuple, trios)) # Re-set them again trios = list(set(trios)) trios = np.array(trios) return trios def fix_biconnections(poss_lost, pairs, gc, pairs_trios): """ Tell me which circle has only 2 connections and what to do about it """ x, y = gc.x, gc.y rmv_trios = [] new_trios = [] # Dictionary of neighbours circs_neighbours = OrderedDict() for i in range(gc.nc): circs_neighbours[i] = [] # Find connections for (i, j) in pairs: circs_neighbours[i].append(j) circs_neighbours[j].append(i) # Set them for i in range(gc.nc): ngh = sorted(circs_neighbours[i]) circs_neighbours[i] = list(set(ngh)) if len(ngh) < 3: # This circle has less than 3 connections # Its connections are: tngh = tuple(ngh) # Their triangles are: alltriosngh = pairs_trios[tngh] # The triangles that dont include it are: notin = [item for item in alltriosngh if i not in item] # Find the triangle where it is contained in for trio in notin: inthisone = gc.triangles[trio].plpol.contains_point((x[i], y[i])) if inthisone: # If it is, this means that this circle is not on the boundary # of the hull, so the issue must be fixed up poss_lost.append(i) # Find its new connection newcon = list(set(trio) - set(tngh))[0] # The two new triangles are: new_trios = [tuple(sorted([i, newcon, item])) for item in tngh] # And this triangle should be removed rmv_trios = [trio] return poss_lost, rmv_trios, new_trios