forked from asher-pembroke/bokeh-tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
filled_contours.py
173 lines (139 loc) · 5.3 KB
/
filled_contours.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
# converts a filled isocontours to a set of polygons for plotting in bokeh
import numpy as np
import matplotlib.pyplot as plt
from webcolors import rgb_to_hex
from bokeh.plotting import figure, show, output_file, hplot
import networkx as nx
from scipy.spatial import Delaunay
def triangulated_graph(polygons):
tri = Delaunay(np.vstack(polygons))
#create ring identifiers and indices into ring vertices
ring_index = []
local_dex = []
for i,poly in enumerate(polygons):
ring_index.append(i+np.zeros(len(poly), dtype=np.int))
local_dex.append(np.arange(len(poly)))
ring_index = np.hstack(ring_index)
local_dex = np.hstack(local_dex)
edges = set()
for simplex in tri.simplices:
edges.add(tuple(sorted([simplex[0], simplex[1]])))
edges.add(tuple(sorted([simplex[0], simplex[2]])))
edges.add(tuple(sorted([simplex[1], simplex[2]])))
#put undirected edges in graph
triangle_graph = nx.Graph()
for e0,e1 in edges:
triangle_graph.add_edge(e0, e1, weight=np.linalg.norm(tri.points[e0]-tri.points[e1]))
# put node data in graph
for i,p in enumerate(tri.points):
triangle_graph.add_node(i,ring=ring_index[i], local = local_dex[i])
return triangle_graph
""" Create a graph of ring islands. some islands will have two bridges joining them """
def create_islands(graph):
#create minimum spanning tree from undirected edges
mst_edges = sorted(list(nx.minimum_spanning_edges(graph,data=True)))
islands = nx.Graph()
for e0, e1, w in mst_edges:
ring0 = graph.node[e0]['ring']
ring1 = graph.node[e1]['ring']
local0, local1 = graph.node[e0]['local'], graph.node[e1]['local']
if ring0 != ring1:
islands.add_edge(ring0, ring1, weight = w,
connection = [e0, e1, local0, local1],
)
return islands
""" Inserts degenerate edge, replacing nodes with new ones
0 -> 1a 1b -> 2 0 -> 1a -> 3b -> 4 -> 5 -> 3a -> 1b -> 2
| |
4 <- 3b 3a <- 5 <- 4
"""
def insert_branch(graph, edge):
for e in edge:
prev, next = graph.predecessors(e), graph.successors(e)
if len(prev) > 0:
graph.add_edge(prev[0],e-.1)
graph.node[e-.1] = graph.node[e]
if len(next) > 0:
graph.add_edge(e+.1,next[0])
graph.node[e+.1] = graph.node[e]
graph.remove_node(e)
# link new nodes. Order won't matter when doing shortest path query
graph.add_edge(edge[0]+.1, edge[1]-.1)
graph.add_edge(edge[0]-.1, edge[1]+.1)
return graph
# create a directed graph from polygons
def polygon_graph(polygons):
g = nx.DiGraph()
shift = 0
for i,poly in enumerate(polygons):
g.add_cycle(range(shift,shift+len(poly)))
shift += len(poly)
return g
"""Get the shortest path that traverses the edges """
def get_merged_path(poly_graph):
#select an edge, remove it and find path that circles back
e1, e0 = poly_graph.edges()[0]
ug = poly_graph.to_undirected()
ug.remove_edge(e0,e1)
# round to get the index of the original vertices
return [int(round(x)) for x in nx.shortest_path(ug,e0)[e1]]
"""Connect each island through one bridge"""
def merge_islands(islands, polygons):
poly_graph = polygon_graph(polygons)
mst_islands = nx.minimum_spanning_tree(islands, 0)
for r_0, r_1 in mst_islands.edges():
e0, e1, local0, local1 = mst_islands[r_0][r_1]['connection']
poly_graph = insert_branch(poly_graph,(e0,e1))
return poly_graph
def filled_contours(p,cn,simplify_threshold = .01):
"""Creates a bokeh plot of filled contours
Args:
p (bokeh.plotting.Figure): Bokeh plot instance
cn (contours): Contours generated from plt.contourf()
simplify_threshold (Optional[float]): Resolution of the output contours in screenspace. Defaults to .01
Returns:
None
"""
for cc in cn.collections:
face_color = np.array(cc.get_facecolor()[0])
color = rgb_to_hex(tuple((255*face_color[:-1]).round().astype(int)))
alpha = face_color[-1]
for path in cc.get_paths():
path.simplify_threshold = simplify_threshold
polygons = path.to_polygons()
if len(polygons) == 1:
p.patch(polygons[0][:,0], polygons[0][:,1], line_alpha = alpha, line_color = color, fill_alpha = alpha, fill_color = color)
else:
vertices = np.vstack(polygons)
graph = triangulated_graph(polygons)
islands = create_islands(graph)
poly_graph = merge_islands(islands,polygons)
merged_path = get_merged_path(poly_graph)
p.patch(vertices[merged_path][:,0],vertices[merged_path][:,1], line_alpha = alpha, line_color = color, fill_alpha = alpha, fill_color = color)
def main(argv):
fig = plt.figure()
ax = fig.add_subplot(121)
N = 200
x = np.linspace(0, 10, N)
y = np.linspace(0, 10, N)
xx, yy = np.meshgrid(x, y)
d = np.sin(xx)*np.cos(2*yy)+.1*yy
d = d + .015*np.random.rand(*xx.shape)
cn = plt.contourf(x,y,d,9, extend = "both")
ax2 = fig.add_subplot(122)
ax2.xlim = [0,10]
ax2.ylim = [0,10]
html_size_without_plots = 1638397 # bytes
html_size_with_p0_only = 1901692 # bytes
html_size_with_p1_only = 2463969 # bytes
output_file('filled_contours.html', title='filled contours')
p0 = figure(x_range=[0, 10], y_range=[0, 10])
filled_contours(p0, cn,.02)
p0.title = str((html_size_with_p0_only-html_size_without_plots)/1000.) + ' Kb'
p1 = figure(x_range=[0, 10], y_range=[0, 10])
p1.image(image = [d], x=[0], y=[0], dw=10, dh=10, palette="Spectral11")
p1.title = str((html_size_with_p1_only-html_size_without_plots)/1000.) + ' Kb'
show(hplot(p0, p1))
if __name__ == '__main__':
import sys
main(sys.argv[1:])