1def iup_segment(coords, rc1, rd1, rc2, rd2): 2 # rc1 = reference coord 1 3 # rd1 = reference delta 1 4 out_arrays = [None, None] 5 for j in 0,1: 6 out_arrays[j] = out = [] 7 x1, x2, d1, d2 = rc1[j], rc2[j], rd1[j], rd2[j] 8 9 10 if x1 == x2: 11 n = len(coords) 12 if d1 == d2: 13 out.extend([d1]*n) 14 else: 15 out.extend([0]*n) 16 continue 17 18 if x1 > x2: 19 x1, x2 = x2, x1 20 d1, d2 = d2, d1 21 22 # x1 < x2 23 scale = (d2 - d1) / (x2 - x1) 24 for pair in coords: 25 x = pair[j] 26 27 if x <= x1: 28 d = d1 29 elif x >= x2: 30 d = d2 31 else: 32 # Interpolate 33 d = d1 + (x - x1) * scale 34 35 out.append(d) 36 37 return zip(*out_arrays) 38 39def iup_contour(delta, coords): 40 assert len(delta) == len(coords) 41 if None not in delta: 42 return delta 43 44 n = len(delta) 45 # indices of points with explicit deltas 46 indices = [i for i,v in enumerate(delta) if v is not None] 47 if not indices: 48 # All deltas are None. Return 0,0 for all. 49 return [(0,0)]*n 50 51 out = [] 52 it = iter(indices) 53 start = next(it) 54 if start != 0: 55 # Initial segment that wraps around 56 i1, i2, ri1, ri2 = 0, start, start, indices[-1] 57 out.extend(iup_segment(coords[i1:i2], coords[ri1], delta[ri1], coords[ri2], delta[ri2])) 58 out.append(delta[start]) 59 for end in it: 60 if end - start > 1: 61 i1, i2, ri1, ri2 = start+1, end, start, end 62 out.extend(iup_segment(coords[i1:i2], coords[ri1], delta[ri1], coords[ri2], delta[ri2])) 63 out.append(delta[end]) 64 start = end 65 if start != n-1: 66 # Final segment that wraps around 67 i1, i2, ri1, ri2 = start+1, n, start, indices[0] 68 out.extend(iup_segment(coords[i1:i2], coords[ri1], delta[ri1], coords[ri2], delta[ri2])) 69 70 assert len(delta) == len(out), (len(delta), len(out)) 71 return out 72 73def iup_delta(delta, coords, ends): 74 assert sorted(ends) == ends and len(coords) == (ends[-1]+1 if ends else 0) + 4 75 n = len(coords) 76 ends = ends + [n-4, n-3, n-2, n-1] 77 out = [] 78 start = 0 79 for end in ends: 80 end += 1 81 contour = iup_contour(delta[start:end], coords[start:end]) 82 out.extend(contour) 83 start = end 84 85 return out 86 87# Optimizer 88 89def can_iup_in_between(deltas, coords, i, j, tolerance): 90 assert j - i >= 2 91 interp = list(iup_segment(coords[i+1:j], coords[i], deltas[i], coords[j], deltas[j])) 92 deltas = deltas[i+1:j] 93 94 assert len(deltas) == len(interp) 95 96 return all(abs(complex(x-p, y-q)) <= tolerance for (x,y),(p,q) in zip(deltas, interp)) 97 98def _iup_contour_bound_forced_set(delta, coords, tolerance=0): 99 """The forced set is a conservative set of points on the contour that must be encoded 100 explicitly (ie. cannot be interpolated). Calculating this set allows for significantly 101 speeding up the dynamic-programming, as well as resolve circularity in DP. 102 103 The set is precise; that is, if an index is in the returned set, then there is no way 104 that IUP can generate delta for that point, given coords and delta. 105 """ 106 assert len(delta) == len(coords) 107 108 forced = set() 109 # Track "last" and "next" points on the contour as we sweep. 110 nd, nc = delta[0], coords[0] 111 ld, lc = delta[-1], coords[-1] 112 for i in range(len(delta)-1, -1, -1): 113 d, c = ld, lc 114 ld, lc = delta[i-1], coords[i-1] 115 116 for j in (0,1): # For X and for Y 117 cj = c[j] 118 dj = d[j] 119 lcj = lc[j] 120 ldj = ld[j] 121 ncj = nc[j] 122 ndj = nd[j] 123 124 if lcj <= ncj: 125 c1, c2 = lcj, ncj 126 d1, d2 = ldj, ndj 127 else: 128 c1, c2 = ncj, lcj 129 d1, d2 = ndj, ldj 130 131 # If coordinate for current point is between coordinate of adjacent 132 # points on the two sides, but the delta for current point is NOT 133 # between delta for those adjacent points (considering tolerance 134 # allowance), then there is no way that current point can be IUP-ed. 135 # Mark it forced. 136 force = False 137 if c1 <= cj <= c2: 138 if not (min(d1,d2)-tolerance <= dj <= max(d1,d2)+tolerance): 139 force = True 140 else: # cj < c1 or c2 < cj 141 if c1 == c2: 142 if d1 == d2: 143 if abs(dj - d1) > tolerance: 144 force = True 145 else: 146 if abs(dj) > tolerance: 147 # Disabled the following because the "d1 == d2" does 148 # check does not take tolerance into consideration... 149 pass # force = True 150 elif d1 != d2: 151 if cj < c1: 152 if dj != d1 and ((dj-tolerance < d1) != (d1 < d2)): 153 force = True 154 else: # c2 < cj 155 if d2 != dj and ((d2 < dj+tolerance) != (d1 < d2)): 156 force = True 157 158 if force: 159 forced.add(i) 160 break 161 162 nd, nc = d, c 163 164 return forced 165 166def _iup_contour_optimize_dp(delta, coords, forced={}, tolerance=0, lookback=None): 167 """Straightforward Dynamic-Programming. For each index i, find least-costly encoding of 168 points 0 to i where i is explicitly encoded. We find this by considering all previous 169 explicit points j and check whether interpolation can fill points between j and i. 170 171 Note that solution always encodes last point explicitly. Higher-level is responsible 172 for removing that restriction. 173 174 As major speedup, we stop looking further whenever we see a "forced" point.""" 175 176 n = len(delta) 177 if lookback is None: 178 lookback = n 179 costs = {-1:0} 180 chain = {-1:None} 181 for i in range(0, n): 182 best_cost = costs[i-1] + 1 183 184 costs[i] = best_cost 185 chain[i] = i - 1 186 187 if i - 1 in forced: 188 continue 189 190 for j in range(i-2, max(i-lookback, -2), -1): 191 192 cost = costs[j] + 1 193 194 if cost < best_cost and can_iup_in_between(delta, coords, j, i, tolerance): 195 costs[i] = best_cost = cost 196 chain[i] = j 197 198 if j in forced: 199 break 200 201 return chain, costs 202 203def _rot_list(l, k): 204 """Rotate list by k items forward. Ie. item at position 0 will be 205 at position k in returned list. Negative k is allowed.""" 206 n = len(l) 207 k %= n 208 if not k: return l 209 return l[n-k:] + l[:n-k] 210 211def _rot_set(s, k, n): 212 k %= n 213 if not k: return s 214 return {(v + k) % n for v in s} 215 216def iup_contour_optimize(delta, coords, tolerance=0.): 217 n = len(delta) 218 219 # Get the easy cases out of the way: 220 221 # If all are within tolerance distance of 0, encode nothing: 222 if all(abs(complex(*p)) <= tolerance for p in delta): 223 return [None] * n 224 225 # If there's exactly one point, return it: 226 if n == 1: 227 return delta 228 229 # If all deltas are exactly the same, return just one (the first one): 230 d0 = delta[0] 231 if all(d0 == d for d in delta): 232 return [d0] + [None] * (n-1) 233 234 # Else, solve the general problem using Dynamic Programming. 235 236 forced = _iup_contour_bound_forced_set(delta, coords, tolerance) 237 # The _iup_contour_optimize_dp() routine returns the optimal encoding 238 # solution given the constraint that the last point is always encoded. 239 # To remove this constraint, we use two different methods, depending on 240 # whether forced set is non-empty or not: 241 242 if forced: 243 # Forced set is non-empty: rotate the contour start point 244 # such that the last point in the list is a forced point. 245 k = (n-1) - max(forced) 246 assert k >= 0 247 248 delta = _rot_list(delta, k) 249 coords = _rot_list(coords, k) 250 forced = _rot_set(forced, k, n) 251 252 chain, costs = _iup_contour_optimize_dp(delta, coords, forced, tolerance) 253 254 # Assemble solution. 255 solution = set() 256 i = n - 1 257 while i is not None: 258 solution.add(i) 259 i = chain[i] 260 assert forced <= solution, (forced, solution) 261 delta = [delta[i] if i in solution else None for i in range(n)] 262 263 delta = _rot_list(delta, -k) 264 else: 265 # Repeat the contour an extra time, solve the 2*n case, then look for solutions of the 266 # circular n-length problem in the solution for 2*n linear case. I cannot prove that 267 # this always produces the optimal solution... 268 chain, costs = _iup_contour_optimize_dp(delta+delta, coords+coords, forced, tolerance, n) 269 best_sol, best_cost = None, n+1 270 271 for start in range(n-1, 2*n-1): 272 # Assemble solution. 273 solution = set() 274 i = start 275 while i > start - n: 276 solution.add(i % n) 277 i = chain[i] 278 if i == start - n: 279 cost = costs[start] - costs[start - n] 280 if cost <= best_cost: 281 best_sol, best_cost = solution, cost 282 283 delta = [delta[i] if i in best_sol else None for i in range(n)] 284 285 286 return delta 287 288def iup_delta_optimize(delta, coords, ends, tolerance=0.): 289 assert sorted(ends) == ends and len(coords) == (ends[-1]+1 if ends else 0) + 4 290 n = len(coords) 291 ends = ends + [n-4, n-3, n-2, n-1] 292 out = [] 293 start = 0 294 for end in ends: 295 contour = iup_contour_optimize(delta[start:end+1], coords[start:end+1], tolerance) 296 assert len(contour) == end - start + 1 297 out.extend(contour) 298 start = end+1 299 300 return out 301