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