1*e1fe3e4aSElliott Hughesfrom fontTools.ttLib.ttGlyphSet import LerpGlyphSet 2*e1fe3e4aSElliott Hughesfrom fontTools.pens.basePen import AbstractPen, BasePen, DecomposingPen 3*e1fe3e4aSElliott Hughesfrom fontTools.pens.pointPen import AbstractPointPen, SegmentToPointPen 4*e1fe3e4aSElliott Hughesfrom fontTools.pens.recordingPen import RecordingPen, DecomposingRecordingPen 5*e1fe3e4aSElliott Hughesfrom fontTools.misc.transform import Transform 6*e1fe3e4aSElliott Hughesfrom collections import defaultdict, deque 7*e1fe3e4aSElliott Hughesfrom math import sqrt, copysign, atan2, pi 8*e1fe3e4aSElliott Hughesfrom enum import Enum 9*e1fe3e4aSElliott Hughesimport itertools 10*e1fe3e4aSElliott Hughes 11*e1fe3e4aSElliott Hughesimport logging 12*e1fe3e4aSElliott Hughes 13*e1fe3e4aSElliott Hugheslog = logging.getLogger("fontTools.varLib.interpolatable") 14*e1fe3e4aSElliott Hughes 15*e1fe3e4aSElliott Hughes 16*e1fe3e4aSElliott Hughesclass InterpolatableProblem: 17*e1fe3e4aSElliott Hughes NOTHING = "nothing" 18*e1fe3e4aSElliott Hughes MISSING = "missing" 19*e1fe3e4aSElliott Hughes OPEN_PATH = "open_path" 20*e1fe3e4aSElliott Hughes PATH_COUNT = "path_count" 21*e1fe3e4aSElliott Hughes NODE_COUNT = "node_count" 22*e1fe3e4aSElliott Hughes NODE_INCOMPATIBILITY = "node_incompatibility" 23*e1fe3e4aSElliott Hughes CONTOUR_ORDER = "contour_order" 24*e1fe3e4aSElliott Hughes WRONG_START_POINT = "wrong_start_point" 25*e1fe3e4aSElliott Hughes KINK = "kink" 26*e1fe3e4aSElliott Hughes UNDERWEIGHT = "underweight" 27*e1fe3e4aSElliott Hughes OVERWEIGHT = "overweight" 28*e1fe3e4aSElliott Hughes 29*e1fe3e4aSElliott Hughes severity = { 30*e1fe3e4aSElliott Hughes MISSING: 1, 31*e1fe3e4aSElliott Hughes OPEN_PATH: 2, 32*e1fe3e4aSElliott Hughes PATH_COUNT: 3, 33*e1fe3e4aSElliott Hughes NODE_COUNT: 4, 34*e1fe3e4aSElliott Hughes NODE_INCOMPATIBILITY: 5, 35*e1fe3e4aSElliott Hughes CONTOUR_ORDER: 6, 36*e1fe3e4aSElliott Hughes WRONG_START_POINT: 7, 37*e1fe3e4aSElliott Hughes KINK: 8, 38*e1fe3e4aSElliott Hughes UNDERWEIGHT: 9, 39*e1fe3e4aSElliott Hughes OVERWEIGHT: 10, 40*e1fe3e4aSElliott Hughes NOTHING: 11, 41*e1fe3e4aSElliott Hughes } 42*e1fe3e4aSElliott Hughes 43*e1fe3e4aSElliott Hughes 44*e1fe3e4aSElliott Hughesdef sort_problems(problems): 45*e1fe3e4aSElliott Hughes """Sort problems by severity, then by glyph name, then by problem message.""" 46*e1fe3e4aSElliott Hughes return dict( 47*e1fe3e4aSElliott Hughes sorted( 48*e1fe3e4aSElliott Hughes problems.items(), 49*e1fe3e4aSElliott Hughes key=lambda _: -min( 50*e1fe3e4aSElliott Hughes ( 51*e1fe3e4aSElliott Hughes (InterpolatableProblem.severity[p["type"]] + p.get("tolerance", 0)) 52*e1fe3e4aSElliott Hughes for p in _[1] 53*e1fe3e4aSElliott Hughes ), 54*e1fe3e4aSElliott Hughes ), 55*e1fe3e4aSElliott Hughes reverse=True, 56*e1fe3e4aSElliott Hughes ) 57*e1fe3e4aSElliott Hughes ) 58*e1fe3e4aSElliott Hughes 59*e1fe3e4aSElliott Hughes 60*e1fe3e4aSElliott Hughesdef rot_list(l, k): 61*e1fe3e4aSElliott Hughes """Rotate list by k items forward. Ie. item at position 0 will be 62*e1fe3e4aSElliott Hughes at position k in returned list. Negative k is allowed.""" 63*e1fe3e4aSElliott Hughes return l[-k:] + l[:-k] 64*e1fe3e4aSElliott Hughes 65*e1fe3e4aSElliott Hughes 66*e1fe3e4aSElliott Hughesclass PerContourPen(BasePen): 67*e1fe3e4aSElliott Hughes def __init__(self, Pen, glyphset=None): 68*e1fe3e4aSElliott Hughes BasePen.__init__(self, glyphset) 69*e1fe3e4aSElliott Hughes self._glyphset = glyphset 70*e1fe3e4aSElliott Hughes self._Pen = Pen 71*e1fe3e4aSElliott Hughes self._pen = None 72*e1fe3e4aSElliott Hughes self.value = [] 73*e1fe3e4aSElliott Hughes 74*e1fe3e4aSElliott Hughes def _moveTo(self, p0): 75*e1fe3e4aSElliott Hughes self._newItem() 76*e1fe3e4aSElliott Hughes self._pen.moveTo(p0) 77*e1fe3e4aSElliott Hughes 78*e1fe3e4aSElliott Hughes def _lineTo(self, p1): 79*e1fe3e4aSElliott Hughes self._pen.lineTo(p1) 80*e1fe3e4aSElliott Hughes 81*e1fe3e4aSElliott Hughes def _qCurveToOne(self, p1, p2): 82*e1fe3e4aSElliott Hughes self._pen.qCurveTo(p1, p2) 83*e1fe3e4aSElliott Hughes 84*e1fe3e4aSElliott Hughes def _curveToOne(self, p1, p2, p3): 85*e1fe3e4aSElliott Hughes self._pen.curveTo(p1, p2, p3) 86*e1fe3e4aSElliott Hughes 87*e1fe3e4aSElliott Hughes def _closePath(self): 88*e1fe3e4aSElliott Hughes self._pen.closePath() 89*e1fe3e4aSElliott Hughes self._pen = None 90*e1fe3e4aSElliott Hughes 91*e1fe3e4aSElliott Hughes def _endPath(self): 92*e1fe3e4aSElliott Hughes self._pen.endPath() 93*e1fe3e4aSElliott Hughes self._pen = None 94*e1fe3e4aSElliott Hughes 95*e1fe3e4aSElliott Hughes def _newItem(self): 96*e1fe3e4aSElliott Hughes self._pen = pen = self._Pen() 97*e1fe3e4aSElliott Hughes self.value.append(pen) 98*e1fe3e4aSElliott Hughes 99*e1fe3e4aSElliott Hughes 100*e1fe3e4aSElliott Hughesclass PerContourOrComponentPen(PerContourPen): 101*e1fe3e4aSElliott Hughes def addComponent(self, glyphName, transformation): 102*e1fe3e4aSElliott Hughes self._newItem() 103*e1fe3e4aSElliott Hughes self.value[-1].addComponent(glyphName, transformation) 104*e1fe3e4aSElliott Hughes 105*e1fe3e4aSElliott Hughes 106*e1fe3e4aSElliott Hughesclass SimpleRecordingPointPen(AbstractPointPen): 107*e1fe3e4aSElliott Hughes def __init__(self): 108*e1fe3e4aSElliott Hughes self.value = [] 109*e1fe3e4aSElliott Hughes 110*e1fe3e4aSElliott Hughes def beginPath(self, identifier=None, **kwargs): 111*e1fe3e4aSElliott Hughes pass 112*e1fe3e4aSElliott Hughes 113*e1fe3e4aSElliott Hughes def endPath(self) -> None: 114*e1fe3e4aSElliott Hughes pass 115*e1fe3e4aSElliott Hughes 116*e1fe3e4aSElliott Hughes def addPoint(self, pt, segmentType=None): 117*e1fe3e4aSElliott Hughes self.value.append((pt, False if segmentType is None else True)) 118*e1fe3e4aSElliott Hughes 119*e1fe3e4aSElliott Hughes 120*e1fe3e4aSElliott Hughesdef vdiff_hypot2(v0, v1): 121*e1fe3e4aSElliott Hughes s = 0 122*e1fe3e4aSElliott Hughes for x0, x1 in zip(v0, v1): 123*e1fe3e4aSElliott Hughes d = x1 - x0 124*e1fe3e4aSElliott Hughes s += d * d 125*e1fe3e4aSElliott Hughes return s 126*e1fe3e4aSElliott Hughes 127*e1fe3e4aSElliott Hughes 128*e1fe3e4aSElliott Hughesdef vdiff_hypot2_complex(v0, v1): 129*e1fe3e4aSElliott Hughes s = 0 130*e1fe3e4aSElliott Hughes for x0, x1 in zip(v0, v1): 131*e1fe3e4aSElliott Hughes d = x1 - x0 132*e1fe3e4aSElliott Hughes s += d.real * d.real + d.imag * d.imag 133*e1fe3e4aSElliott Hughes # This does the same but seems to be slower: 134*e1fe3e4aSElliott Hughes # s += (d * d.conjugate()).real 135*e1fe3e4aSElliott Hughes return s 136*e1fe3e4aSElliott Hughes 137*e1fe3e4aSElliott Hughes 138*e1fe3e4aSElliott Hughesdef matching_cost(G, matching): 139*e1fe3e4aSElliott Hughes return sum(G[i][j] for i, j in enumerate(matching)) 140*e1fe3e4aSElliott Hughes 141*e1fe3e4aSElliott Hughes 142*e1fe3e4aSElliott Hughesdef min_cost_perfect_bipartite_matching_scipy(G): 143*e1fe3e4aSElliott Hughes n = len(G) 144*e1fe3e4aSElliott Hughes rows, cols = linear_sum_assignment(G) 145*e1fe3e4aSElliott Hughes assert (rows == list(range(n))).all() 146*e1fe3e4aSElliott Hughes return list(cols), matching_cost(G, cols) 147*e1fe3e4aSElliott Hughes 148*e1fe3e4aSElliott Hughes 149*e1fe3e4aSElliott Hughesdef min_cost_perfect_bipartite_matching_munkres(G): 150*e1fe3e4aSElliott Hughes n = len(G) 151*e1fe3e4aSElliott Hughes cols = [None] * n 152*e1fe3e4aSElliott Hughes for row, col in Munkres().compute(G): 153*e1fe3e4aSElliott Hughes cols[row] = col 154*e1fe3e4aSElliott Hughes return cols, matching_cost(G, cols) 155*e1fe3e4aSElliott Hughes 156*e1fe3e4aSElliott Hughes 157*e1fe3e4aSElliott Hughesdef min_cost_perfect_bipartite_matching_bruteforce(G): 158*e1fe3e4aSElliott Hughes n = len(G) 159*e1fe3e4aSElliott Hughes 160*e1fe3e4aSElliott Hughes if n > 6: 161*e1fe3e4aSElliott Hughes raise Exception("Install Python module 'munkres' or 'scipy >= 0.17.0'") 162*e1fe3e4aSElliott Hughes 163*e1fe3e4aSElliott Hughes # Otherwise just brute-force 164*e1fe3e4aSElliott Hughes permutations = itertools.permutations(range(n)) 165*e1fe3e4aSElliott Hughes best = list(next(permutations)) 166*e1fe3e4aSElliott Hughes best_cost = matching_cost(G, best) 167*e1fe3e4aSElliott Hughes for p in permutations: 168*e1fe3e4aSElliott Hughes cost = matching_cost(G, p) 169*e1fe3e4aSElliott Hughes if cost < best_cost: 170*e1fe3e4aSElliott Hughes best, best_cost = list(p), cost 171*e1fe3e4aSElliott Hughes return best, best_cost 172*e1fe3e4aSElliott Hughes 173*e1fe3e4aSElliott Hughes 174*e1fe3e4aSElliott Hughestry: 175*e1fe3e4aSElliott Hughes from scipy.optimize import linear_sum_assignment 176*e1fe3e4aSElliott Hughes 177*e1fe3e4aSElliott Hughes min_cost_perfect_bipartite_matching = min_cost_perfect_bipartite_matching_scipy 178*e1fe3e4aSElliott Hughesexcept ImportError: 179*e1fe3e4aSElliott Hughes try: 180*e1fe3e4aSElliott Hughes from munkres import Munkres 181*e1fe3e4aSElliott Hughes 182*e1fe3e4aSElliott Hughes min_cost_perfect_bipartite_matching = ( 183*e1fe3e4aSElliott Hughes min_cost_perfect_bipartite_matching_munkres 184*e1fe3e4aSElliott Hughes ) 185*e1fe3e4aSElliott Hughes except ImportError: 186*e1fe3e4aSElliott Hughes min_cost_perfect_bipartite_matching = ( 187*e1fe3e4aSElliott Hughes min_cost_perfect_bipartite_matching_bruteforce 188*e1fe3e4aSElliott Hughes ) 189*e1fe3e4aSElliott Hughes 190*e1fe3e4aSElliott Hughes 191*e1fe3e4aSElliott Hughesdef contour_vector_from_stats(stats): 192*e1fe3e4aSElliott Hughes # Don't change the order of items here. 193*e1fe3e4aSElliott Hughes # It's okay to add to the end, but otherwise, other 194*e1fe3e4aSElliott Hughes # code depends on it. Search for "covariance". 195*e1fe3e4aSElliott Hughes size = sqrt(abs(stats.area)) 196*e1fe3e4aSElliott Hughes return ( 197*e1fe3e4aSElliott Hughes copysign((size), stats.area), 198*e1fe3e4aSElliott Hughes stats.meanX, 199*e1fe3e4aSElliott Hughes stats.meanY, 200*e1fe3e4aSElliott Hughes stats.stddevX * 2, 201*e1fe3e4aSElliott Hughes stats.stddevY * 2, 202*e1fe3e4aSElliott Hughes stats.correlation * size, 203*e1fe3e4aSElliott Hughes ) 204*e1fe3e4aSElliott Hughes 205*e1fe3e4aSElliott Hughes 206*e1fe3e4aSElliott Hughesdef matching_for_vectors(m0, m1): 207*e1fe3e4aSElliott Hughes n = len(m0) 208*e1fe3e4aSElliott Hughes 209*e1fe3e4aSElliott Hughes identity_matching = list(range(n)) 210*e1fe3e4aSElliott Hughes 211*e1fe3e4aSElliott Hughes costs = [[vdiff_hypot2(v0, v1) for v1 in m1] for v0 in m0] 212*e1fe3e4aSElliott Hughes ( 213*e1fe3e4aSElliott Hughes matching, 214*e1fe3e4aSElliott Hughes matching_cost, 215*e1fe3e4aSElliott Hughes ) = min_cost_perfect_bipartite_matching(costs) 216*e1fe3e4aSElliott Hughes identity_cost = sum(costs[i][i] for i in range(n)) 217*e1fe3e4aSElliott Hughes return matching, matching_cost, identity_cost 218*e1fe3e4aSElliott Hughes 219*e1fe3e4aSElliott Hughes 220*e1fe3e4aSElliott Hughesdef points_characteristic_bits(points): 221*e1fe3e4aSElliott Hughes bits = 0 222*e1fe3e4aSElliott Hughes for pt, b in reversed(points): 223*e1fe3e4aSElliott Hughes bits = (bits << 1) | b 224*e1fe3e4aSElliott Hughes return bits 225*e1fe3e4aSElliott Hughes 226*e1fe3e4aSElliott Hughes 227*e1fe3e4aSElliott Hughes_NUM_ITEMS_PER_POINTS_COMPLEX_VECTOR = 4 228*e1fe3e4aSElliott Hughes 229*e1fe3e4aSElliott Hughes 230*e1fe3e4aSElliott Hughesdef points_complex_vector(points): 231*e1fe3e4aSElliott Hughes vector = [] 232*e1fe3e4aSElliott Hughes if not points: 233*e1fe3e4aSElliott Hughes return vector 234*e1fe3e4aSElliott Hughes points = [complex(*pt) for pt, _ in points] 235*e1fe3e4aSElliott Hughes n = len(points) 236*e1fe3e4aSElliott Hughes assert _NUM_ITEMS_PER_POINTS_COMPLEX_VECTOR == 4 237*e1fe3e4aSElliott Hughes points.extend(points[: _NUM_ITEMS_PER_POINTS_COMPLEX_VECTOR - 1]) 238*e1fe3e4aSElliott Hughes while len(points) < _NUM_ITEMS_PER_POINTS_COMPLEX_VECTOR: 239*e1fe3e4aSElliott Hughes points.extend(points[: _NUM_ITEMS_PER_POINTS_COMPLEX_VECTOR - 1]) 240*e1fe3e4aSElliott Hughes for i in range(n): 241*e1fe3e4aSElliott Hughes # The weights are magic numbers. 242*e1fe3e4aSElliott Hughes 243*e1fe3e4aSElliott Hughes # The point itself 244*e1fe3e4aSElliott Hughes p0 = points[i] 245*e1fe3e4aSElliott Hughes vector.append(p0) 246*e1fe3e4aSElliott Hughes 247*e1fe3e4aSElliott Hughes # The vector to the next point 248*e1fe3e4aSElliott Hughes p1 = points[i + 1] 249*e1fe3e4aSElliott Hughes d0 = p1 - p0 250*e1fe3e4aSElliott Hughes vector.append(d0 * 3) 251*e1fe3e4aSElliott Hughes 252*e1fe3e4aSElliott Hughes # The turn vector 253*e1fe3e4aSElliott Hughes p2 = points[i + 2] 254*e1fe3e4aSElliott Hughes d1 = p2 - p1 255*e1fe3e4aSElliott Hughes vector.append(d1 - d0) 256*e1fe3e4aSElliott Hughes 257*e1fe3e4aSElliott Hughes # The angle to the next point, as a cross product; 258*e1fe3e4aSElliott Hughes # Square root of, to match dimentionality of distance. 259*e1fe3e4aSElliott Hughes cross = d0.real * d1.imag - d0.imag * d1.real 260*e1fe3e4aSElliott Hughes cross = copysign(sqrt(abs(cross)), cross) 261*e1fe3e4aSElliott Hughes vector.append(cross * 4) 262*e1fe3e4aSElliott Hughes 263*e1fe3e4aSElliott Hughes return vector 264*e1fe3e4aSElliott Hughes 265*e1fe3e4aSElliott Hughes 266*e1fe3e4aSElliott Hughesdef add_isomorphisms(points, isomorphisms, reverse): 267*e1fe3e4aSElliott Hughes reference_bits = points_characteristic_bits(points) 268*e1fe3e4aSElliott Hughes n = len(points) 269*e1fe3e4aSElliott Hughes 270*e1fe3e4aSElliott Hughes # if points[0][0] == points[-1][0]: 271*e1fe3e4aSElliott Hughes # abort 272*e1fe3e4aSElliott Hughes 273*e1fe3e4aSElliott Hughes if reverse: 274*e1fe3e4aSElliott Hughes points = points[::-1] 275*e1fe3e4aSElliott Hughes bits = points_characteristic_bits(points) 276*e1fe3e4aSElliott Hughes else: 277*e1fe3e4aSElliott Hughes bits = reference_bits 278*e1fe3e4aSElliott Hughes 279*e1fe3e4aSElliott Hughes vector = points_complex_vector(points) 280*e1fe3e4aSElliott Hughes 281*e1fe3e4aSElliott Hughes assert len(vector) % n == 0 282*e1fe3e4aSElliott Hughes mult = len(vector) // n 283*e1fe3e4aSElliott Hughes mask = (1 << n) - 1 284*e1fe3e4aSElliott Hughes 285*e1fe3e4aSElliott Hughes for i in range(n): 286*e1fe3e4aSElliott Hughes b = ((bits << (n - i)) & mask) | (bits >> i) 287*e1fe3e4aSElliott Hughes if b == reference_bits: 288*e1fe3e4aSElliott Hughes isomorphisms.append( 289*e1fe3e4aSElliott Hughes (rot_list(vector, -i * mult), n - 1 - i if reverse else i, reverse) 290*e1fe3e4aSElliott Hughes ) 291*e1fe3e4aSElliott Hughes 292*e1fe3e4aSElliott Hughes 293*e1fe3e4aSElliott Hughesdef find_parents_and_order(glyphsets, locations): 294*e1fe3e4aSElliott Hughes parents = [None] + list(range(len(glyphsets) - 1)) 295*e1fe3e4aSElliott Hughes order = list(range(len(glyphsets))) 296*e1fe3e4aSElliott Hughes if locations: 297*e1fe3e4aSElliott Hughes # Order base master first 298*e1fe3e4aSElliott Hughes bases = (i for i, l in enumerate(locations) if all(v == 0 for v in l.values())) 299*e1fe3e4aSElliott Hughes if bases: 300*e1fe3e4aSElliott Hughes base = next(bases) 301*e1fe3e4aSElliott Hughes logging.info("Base master index %s, location %s", base, locations[base]) 302*e1fe3e4aSElliott Hughes else: 303*e1fe3e4aSElliott Hughes base = 0 304*e1fe3e4aSElliott Hughes logging.warning("No base master location found") 305*e1fe3e4aSElliott Hughes 306*e1fe3e4aSElliott Hughes # Form a minimum spanning tree of the locations 307*e1fe3e4aSElliott Hughes try: 308*e1fe3e4aSElliott Hughes from scipy.sparse.csgraph import minimum_spanning_tree 309*e1fe3e4aSElliott Hughes 310*e1fe3e4aSElliott Hughes graph = [[0] * len(locations) for _ in range(len(locations))] 311*e1fe3e4aSElliott Hughes axes = set() 312*e1fe3e4aSElliott Hughes for l in locations: 313*e1fe3e4aSElliott Hughes axes.update(l.keys()) 314*e1fe3e4aSElliott Hughes axes = sorted(axes) 315*e1fe3e4aSElliott Hughes vectors = [tuple(l.get(k, 0) for k in axes) for l in locations] 316*e1fe3e4aSElliott Hughes for i, j in itertools.combinations(range(len(locations)), 2): 317*e1fe3e4aSElliott Hughes graph[i][j] = vdiff_hypot2(vectors[i], vectors[j]) 318*e1fe3e4aSElliott Hughes 319*e1fe3e4aSElliott Hughes tree = minimum_spanning_tree(graph) 320*e1fe3e4aSElliott Hughes rows, cols = tree.nonzero() 321*e1fe3e4aSElliott Hughes graph = defaultdict(set) 322*e1fe3e4aSElliott Hughes for row, col in zip(rows, cols): 323*e1fe3e4aSElliott Hughes graph[row].add(col) 324*e1fe3e4aSElliott Hughes graph[col].add(row) 325*e1fe3e4aSElliott Hughes 326*e1fe3e4aSElliott Hughes # Traverse graph from the base and assign parents 327*e1fe3e4aSElliott Hughes parents = [None] * len(locations) 328*e1fe3e4aSElliott Hughes order = [] 329*e1fe3e4aSElliott Hughes visited = set() 330*e1fe3e4aSElliott Hughes queue = deque([base]) 331*e1fe3e4aSElliott Hughes while queue: 332*e1fe3e4aSElliott Hughes i = queue.popleft() 333*e1fe3e4aSElliott Hughes visited.add(i) 334*e1fe3e4aSElliott Hughes order.append(i) 335*e1fe3e4aSElliott Hughes for j in sorted(graph[i]): 336*e1fe3e4aSElliott Hughes if j not in visited: 337*e1fe3e4aSElliott Hughes parents[j] = i 338*e1fe3e4aSElliott Hughes queue.append(j) 339*e1fe3e4aSElliott Hughes 340*e1fe3e4aSElliott Hughes except ImportError: 341*e1fe3e4aSElliott Hughes pass 342*e1fe3e4aSElliott Hughes 343*e1fe3e4aSElliott Hughes log.info("Parents: %s", parents) 344*e1fe3e4aSElliott Hughes log.info("Order: %s", order) 345*e1fe3e4aSElliott Hughes return parents, order 346*e1fe3e4aSElliott Hughes 347*e1fe3e4aSElliott Hughes 348*e1fe3e4aSElliott Hughesdef transform_from_stats(stats, inverse=False): 349*e1fe3e4aSElliott Hughes # https://cookierobotics.com/007/ 350*e1fe3e4aSElliott Hughes a = stats.varianceX 351*e1fe3e4aSElliott Hughes b = stats.covariance 352*e1fe3e4aSElliott Hughes c = stats.varianceY 353*e1fe3e4aSElliott Hughes 354*e1fe3e4aSElliott Hughes delta = (((a - c) * 0.5) ** 2 + b * b) ** 0.5 355*e1fe3e4aSElliott Hughes lambda1 = (a + c) * 0.5 + delta # Major eigenvalue 356*e1fe3e4aSElliott Hughes lambda2 = (a + c) * 0.5 - delta # Minor eigenvalue 357*e1fe3e4aSElliott Hughes theta = atan2(lambda1 - a, b) if b != 0 else (pi * 0.5 if a < c else 0) 358*e1fe3e4aSElliott Hughes trans = Transform() 359*e1fe3e4aSElliott Hughes 360*e1fe3e4aSElliott Hughes if lambda2 < 0: 361*e1fe3e4aSElliott Hughes # XXX This is a hack. 362*e1fe3e4aSElliott Hughes # The problem is that the covariance matrix is singular. 363*e1fe3e4aSElliott Hughes # This happens when the contour is a line, or a circle. 364*e1fe3e4aSElliott Hughes # In that case, the covariance matrix is not a good 365*e1fe3e4aSElliott Hughes # representation of the contour. 366*e1fe3e4aSElliott Hughes # We should probably detect this earlier and avoid 367*e1fe3e4aSElliott Hughes # computing the covariance matrix in the first place. 368*e1fe3e4aSElliott Hughes # But for now, we just avoid the division by zero. 369*e1fe3e4aSElliott Hughes lambda2 = 0 370*e1fe3e4aSElliott Hughes 371*e1fe3e4aSElliott Hughes if inverse: 372*e1fe3e4aSElliott Hughes trans = trans.translate(-stats.meanX, -stats.meanY) 373*e1fe3e4aSElliott Hughes trans = trans.rotate(-theta) 374*e1fe3e4aSElliott Hughes trans = trans.scale(1 / sqrt(lambda1), 1 / sqrt(lambda2)) 375*e1fe3e4aSElliott Hughes else: 376*e1fe3e4aSElliott Hughes trans = trans.scale(sqrt(lambda1), sqrt(lambda2)) 377*e1fe3e4aSElliott Hughes trans = trans.rotate(theta) 378*e1fe3e4aSElliott Hughes trans = trans.translate(stats.meanX, stats.meanY) 379*e1fe3e4aSElliott Hughes 380*e1fe3e4aSElliott Hughes return trans 381