# Ticket #11564: polyhedron_unfold.py

File polyhedron_unfold.py, 13.8 KB (added by , 8 years ago) |
---|

Line | |
---|---|

1 | from sage.geometry.polyhedra import cyclic_sort_vertices_2d |

2 | from sage.plot.line import line |

3 | from sage.plot.polygon import polygon |

4 | from math import sqrt |

5 | from copy import copy |

6 | from random import random |

7 | |

8 | ########################################################## |

9 | # |

10 | # Returns the dot product between two vectors |

11 | # |

12 | ########################################################## |

13 | |

14 | def dot_prod(v1,v2): |

15 | dot=0.0 |

16 | for i in range(len(v1)): |

17 | dot+=v1[i]*v2[i] |

18 | return float(dot) |

19 | |

20 | ########################################################## |

21 | # |

22 | # Returns the length of a vector |

23 | # |

24 | ########################################################## |

25 | |

26 | def length(v): |

27 | d=0.0 |

28 | for i in v: |

29 | d += i**2.0 |

30 | return float(sqrt(d)) |

31 | |

32 | ########################################################## |

33 | # |

34 | # Returns the length of an edge |

35 | # |

36 | ########################################################## |

37 | |

38 | def edge_length(edge): |

39 | d=0.0 |

40 | for i in range(len(edge[0])): |

41 | d += (edge[1][i]-edge[0][i])**2.0 |

42 | return float(sqrt(d)) |

43 | |

44 | ########################################################## |

45 | # |

46 | # Returns a vector between two points |

47 | # |

48 | ########################################################## |

49 | |

50 | def vector(p1,p2): |

51 | v=[] |

52 | for i in range(len(p1)): |

53 | v.append(p2[i]-p1[i]) |

54 | return v |

55 | |

56 | ########################################################## |

57 | # |

58 | # Returns a vector between two points of an edge |

59 | # |

60 | ########################################################## |

61 | |

62 | def edge_vector(edge): |

63 | v=[] |

64 | return vector(edge[0],edge[1]) |

65 | |

66 | ########################################################## |

67 | # |

68 | # Returns the cosine of the angle between two vectors |

69 | # |

70 | ########################################################## |

71 | |

72 | def cos_angle(v1,v2): |

73 | dot = dot_prod(v1,v2) |

74 | d1 = length(v1) |

75 | d2 = length(v2) |

76 | return float(dot/(d1*d2)) |

77 | |

78 | ########################################################## |

79 | # |

80 | # Returns a map of the face from 3D to 2D |

81 | # |

82 | ########################################################## |

83 | |

84 | def face_to_plane(face_edges,face_edge_index=None,map_edge=None,alt=1): |

85 | """ |

86 | Returns a map of the 3D face to 2D. For use in poly_unfold. |

87 | |

88 | INPUT: |

89 | |

90 | - ``face_edges`` -- A list of edges which define a face in 3D. |

91 | - ``face_edge_index`` -- An integer which is the index of the edge in face_edges which is being attached. |

92 | - ``map_edge`` -- A list of two anchor points where the edge at face_edge_index will be attached. |

93 | - ``alt`` -- Either 1 or -1. A face can be mapped in two different ways given two anchor points. This picks one. |

94 | |

95 | OUTPUT: |

96 | |

97 | A list of edges which define the mapping of the face to 2D. |

98 | |

99 | EXAMPLES:: |

100 | |

101 | sage: face_edges = [[(-1,-1,1),(-1,1,1)],[(-1,1,1),(1,1,1)],[(1,1,1),(1,-1,1)],[(1,-1,1),(-1,-1,1)]] |

102 | sage: map_edge = [(0,0),(1,0)] |

103 | sage: face_to_plane(face_edges,0,map_edge) |

104 | [[(0, 0), (1, 0)], [(2.0, 0.0), (2.0, 2.0)], [(2.0, 2.0), (0.0, 2.0)], [(0.0, 2.0), (0, 0)]] |

105 | """ |

106 | n = len(face_edges) |

107 | |

108 | face_map = n*[None] |

109 | |

110 | if face_edge_index is None: |

111 | face_edge_index = 0 |

112 | |

113 | face_edge = face_edges[face_edge_index] |

114 | |

115 | if map_edge is None: |

116 | d = edge_length(face_edge) |

117 | map_edge = [(0,0),(d,0)] |

118 | |

119 | face_map[face_edge_index] = map_edge |

120 | |

121 | for i in range(face_edge_index-1,face_edge_index-n,-1): |

122 | |

123 | unit_dir = edge_vector(face_map[i+1]) |

124 | d = edge_length(face_map[i+1]) |

125 | for j in range(2): |

126 | unit_dir[j]/=d |

127 | |

128 | v1 = edge_vector(face_edges[i+1]) |

129 | v2 = edge_vector(face_edges[i]) |

130 | |

131 | d = length(v2) |

132 | |

133 | c=cos_angle(v1,v2) |

134 | s=alt*sqrt(1.0-c*c) |

135 | |

136 | x = (-c*unit_dir[0] - s*unit_dir[1])*d + face_map[i+1][0][0] |

137 | y = (s*unit_dir[0] - c*unit_dir[1])*d + face_map[i+1][0][1] |

138 | |

139 | face_map[i] = [(x,y),face_map[i+1][0]] |

140 | |

141 | return face_map |

142 | |

143 | ########################################################## |

144 | # |

145 | # Returns face adjacency data and a list of face edges |

146 | # |

147 | ########################################################## |

148 | |

149 | def face_adj_data(poly): |

150 | """ |

151 | Returns a dictionary of dictionaries containing information about face adjacencies |

152 | and a list of faces (a face here is a list of edges). |

153 | |

154 | INPUT: |

155 | |

156 | - ``poly`` -- A 3D sage Polyhedron. |

157 | |

158 | OUTPUT: |

159 | |

160 | A dictionary of dictionaries and a list of faces. |

161 | The dicitionary of dictionaries has as keys the index of the face (in faces, call it k1), and as |

162 | value a dictionary with the index of the adjacent face (call it k2) as key and the index of the edge |

163 | (in the list of edges for the face at index k1) as value. |

164 | |

165 | |

166 | EXAMPLES:: |

167 | |

168 | sage: cube=polytopes.n_cube(3) |

169 | sage: data,faces = face_adj_data(cube) |

170 | sage: data |

171 | {0: {1: 2, 2: 3, 3: 1, 5: 0}, 1: {0: 2, 2: 3, 3: 1, 4: 0}, 2: {0: 2, 1: 3, 4: 0, 5: 1}, |

172 | 3: {0: 2, 1: 3, 4: 0, 5: 1}, 4: {1: 2, 2: 3, 3: 1, 5: 0}, 5: {0: 2, 2: 3, 3: 1, 4: 0}} |

173 | """ |

174 | faces = [] |

175 | face_lattice_ls = poly.face_lattice().level_sets() |

176 | n = len(face_lattice_ls) |

177 | |

178 | for i in range(3,n-1): |

179 | faces += face_lattice_ls[i] |

180 | |

181 | data = {} |

182 | n = len(faces) |

183 | |

184 | for i in range(n): |

185 | data.update({i:{}}) |

186 | f = [k for k in faces[i].element.Vrepresentation()] |

187 | f = cyclic_sort_vertices_2d(f) |

188 | faces[i]=[] |

189 | for j in range(-1,len(f)-1): |

190 | faces[i].append([f[j],f[j+1]]) |

191 | |

192 | for i in range(n): |

193 | for j in range(i+1,n): |

194 | sfi = map(set,faces[i]) |

195 | sfj = map(set,faces[j]) |

196 | for a in sfi: |

197 | for b in sfj: |

198 | if a == b: |

199 | data[i].update({j:sfi.index(a)}) |

200 | data[j].update({i:sfj.index(b)}) |

201 | break |

202 | |

203 | return data,faces |

204 | |

205 | ############################################################## |

206 | # |

207 | # Returns the coordinates of the centroid of a list of edges |

208 | # |

209 | ############################################################## |

210 | |

211 | def centroid(shape): |

212 | """ |

213 | Returns the centroid of a list of edges. |

214 | |

215 | INPUT: |

216 | |

217 | - ``shape`` -- A list of edges. |

218 | |

219 | OUTPUT: |

220 | |

221 | A point which is the centroid of shape. |

222 | |

223 | EXAMPLES:: |

224 | |

225 | sage: shape = [[[0,0],[1,0]],[[1,0],[1,1]],[[1,1],[0,1]],[[0,1],[0,0]]] |

226 | sage: centroid(shape) |

227 | [0.5, 0.5] |

228 | """ |

229 | n = len(shape[0][0]) |

230 | c=n*[0.0] |

231 | for edge in shape: |

232 | for i in edge: |

233 | for j in range(n): |

234 | c[j] += i[j] |

235 | for i in range(n): |

236 | c[i] /= (2.0*len(shape)) |

237 | return c |

238 | |

239 | ########################################################## |

240 | # |

241 | # Function to determine if point ``p`` is in net ``net`` |

242 | # |

243 | ########################################################## |

244 | |

245 | def point_in_net(p,net): |

246 | """ |

247 | Determines wether or not the point is in the net. Works in 2D |

248 | |

249 | INPUT: |

250 | |

251 | - ``p`` -- A 2D point. |

252 | - ``net`` -- A net, which is a list of 2D edges. |

253 | |

254 | OUTPUT: |

255 | |

256 | True or False, depending on wether the point is in or out of the net. |

257 | |

258 | EXAMPLES:: |

259 | |

260 | sage: shape = [[[0,0],[1,0]],[[1,0],[1,1]],[[1,1],[0,1]],[[0,1],[0,0]]] |

261 | sage: point_in_net(centroid(shape),shape) |

262 | True |

263 | """ |

264 | eps=0.01 |

265 | x = p[0] |

266 | y = p[1] |

267 | count = 0 |

268 | for edge in net: |

269 | if (y-edge[0][1])*(y-edge[1][1]) < 0: |

270 | t = (y - edge[0][1])/(edge[1][1]-edge[0][1]) |

271 | xval = edge[0][0] + (edge[1][0]-edge[0][0])*t |

272 | if x < xval: |

273 | count += 1 |

274 | if count%2 == 1: |

275 | return True |

276 | return False |

277 | |

278 | ################################################################### |

279 | # |

280 | # Function to determine if edge ``e1`` is the same as edge ``e2`` |

281 | # |

282 | ################################################################### |

283 | |

284 | def edge_equal(e1,e2): |

285 | eps = 0.01 |

286 | if abs(e1[0][0]-e2[0][0]) < eps and abs(e1[0][1]-e2[0][1]) < eps: |

287 | if abs(e1[1][0]-e2[1][0]) < eps and abs(e1[1][1]-e2[1][1]) < eps: |

288 | return True |

289 | if abs(e1[1][0]-e2[0][0]) < eps and abs(e1[1][1]-e2[0][1]) < eps: |

290 | if abs(e1[0][0]-e2[1][0]) < eps and abs(e1[0][1]-e2[1][1]) < eps: |

291 | return True |

292 | return False |

293 | |

294 | ########################################################## |

295 | # |

296 | # Function to determine if map ``m`` is in net ``net`` |

297 | # |

298 | ########################################################## |

299 | |

300 | def map_intersects_net(m,net): |

301 | """ |

302 | Determines wether or not the map intersects the net. Works in 2D. |

303 | |

304 | INPUT: |

305 | |

306 | - ``m`` -- A map (a list of edges). |

307 | - ``net`` -- A net (a list of edges). |

308 | |

309 | OUTPUT: |

310 | |

311 | True or False, depending on wether the map intersects the net or not. |

312 | |

313 | EXAMPLES:: |

314 | |

315 | sage: net = [[0,0],[1,0]],[[1,0],[1,1]],[[1,1],[0,1]],[[0,1],[0,0]]] |

316 | sage: m = [[[-0.5,-0.5],[0.5,-0.5]],[[0.5,-0.5],[0.5,0.5]],[[0.5,0.5],[-0.5,0.5]],[[-0.5,0.5],[-0.5,-0.5]]] |

317 | sage: map_intersects_net(m,net) |

318 | True |

319 | """ |

320 | if point_in_net(centroid(m),net): |

321 | return True |

322 | |

323 | eps = 0.01 |

324 | count = 0 |

325 | for j in m: |

326 | x2=j[1][0]-j[0][0] |

327 | y2=j[1][1]-j[0][1] |

328 | |

329 | for i in net: |

330 | |

331 | if edge_equal(i,j): |

332 | count += 1 |

333 | if count == 2: |

334 | return True |

335 | |

336 | x1 = i[1][0]-i[0][0] |

337 | y1 = i[1][1]-i[0][1] |

338 | |

339 | a = (-y1*x2 + x1*y2) |

340 | |

341 | if abs(a) <= eps: |

342 | continue |

343 | |

344 | x3 = i[0][0]-j[0][0] |

345 | y3 = i[0][1]-j[0][1] |

346 | |

347 | a = (-y1*x3 + x1*y3)/a |

348 | b = -(-y2*x3 + x2*y3)/(-y2*x1 + x2*y1) |

349 | |

350 | if (a > eps and a < 1-eps) and (b > -eps and b < 1+eps): |

351 | return True |

352 | elif (b > eps and b < 1-eps) and (a > -eps and a < 1+eps): |

353 | return True |

354 | |

355 | return False |

356 | |

357 | ########################################################## |

358 | # |

359 | # Function to unfold a polyhedron |

360 | # |

361 | ########################################################## |

362 | |

363 | def unfold_poly(data,faces,x=0): |

364 | """ |

365 | Attempts to unfold the polyhedron defined by the input. |

366 | |

367 | INPUT: |

368 | |

369 | - ``data`` -- A dictionary of dictionaries generated by face_adj_data. |

370 | - ``faces`` -- A list of edges generated by face_adj_data. |

371 | - ``x`` -- The index of the starting face. |

372 | |

373 | OUTPUT: |

374 | |

375 | A plot of the unfolding if possible. |

376 | |

377 | """ |

378 | n = len(data) |

379 | #print 'Attempt',x+1 |

380 | |

381 | if x >= n: |

382 | print 'Unable to unfold polyhedron' |

383 | return |

384 | |

385 | previous_faces = [] |

386 | previous_faces.append(x) |

387 | |

388 | previous_maps = n*[None] |

389 | previous_maps[x] = face_to_plane(faces[x]) |

390 | |

391 | neighbour_data = {x:copy(data[x])} |

392 | |

393 | net = copy(previous_maps[x]) |

394 | |

395 | altrn = n*[None] |

396 | altrn[x] = -1 |

397 | |

398 | while len(previous_faces) < n: |

399 | |

400 | failed = True |

401 | |

402 | for old in neighbour_data.keys(): |

403 | |

404 | for new in neighbour_data[old].keys(): |

405 | |

406 | a = data[old][new] |

407 | b = data[new][old] |

408 | |

409 | alt = altrn[old] |

410 | |

411 | old_map = previous_maps[old] |

412 | new_map = face_to_plane(faces[new],b,old_map[a],alt) |

413 | |

414 | if not map_intersects_net(new_map,net): |

415 | |

416 | failed = False |

417 | previous_maps[new] = new_map |

418 | altrn[new] = -alt |

419 | new_data = {new:copy(data[new])} |

420 | previous_faces.append(new) |

421 | |

422 | for i in neighbour_data: |

423 | if new in neighbour_data[i]: |

424 | neighbour_data[i].pop(new) |

425 | if i in new_data[new]: |

426 | new_data[new].pop(i) |

427 | neighbour_data.update(new_data) |

428 | |

429 | net += new_map |

430 | |

431 | if failed: |

432 | unfold_poly(data,faces,x+1) |

433 | return |

434 | |

435 | return plot_maps(previous_maps) |

436 | |

437 | ########################################################## |

438 | # |

439 | # Driver function which calls unfold_poly |

440 | # |

441 | ########################################################## |

442 | |

443 | def UnfoldPolyhedron(poly): |

444 | """ |

445 | Driver function which calls unfold_poly. Attempts to plot the polyhedron's |

446 | unfolding (i.e. net). |

447 | |

448 | INPUT: |

449 | |

450 | - ``poly`` -- A 3D sage Polyhedron (type <class 'sage.geometry.polyhedra.Polyhedron'>). |

451 | |

452 | OUTPUT: |

453 | |

454 | A plot of the unfolding if possible. |

455 | |

456 | EXAMPLES:: |

457 | |

458 | sage: gr = polytopes.great_rhombicuboctahedron() |

459 | sage: UnfoldPolyhedron(gr) |

460 | sage: |

461 | |

462 | """ |

463 | data, faces = face_adj_data(poly) |

464 | #print len(faces),'faces' |

465 | return unfold_poly(data,faces) |

466 | |

467 | ########################################################## |

468 | # |

469 | # Plots the list of 2D face maps |

470 | # |

471 | ########################################################## |

472 | |

473 | def plot_maps(maps): |

474 | """ |

475 | Plots a list of 2D mappings using polygons and lines. |

476 | |

477 | INPUT: |

478 | |

479 | - ``maps`` -- A list of faces, which are lists of edges. |

480 | |

481 | OUTPUT: |

482 | |

483 | A plot of filled polygons and lines. |

484 | |

485 | EXAMPLES:: |

486 | |

487 | sage: maps = [[[[0, 0], [1, 0]], [[1, 0], [1, 1]], [[1, 1], [0, 1]], [[0, 1], [0, 0]]]] |

488 | sage: plot_maps(maps) |

489 | sage: |

490 | |

491 | """ |

492 | f=None |

493 | g=None |

494 | c1=(random(),random(),random()) |

495 | c2=(random(),random(),random()) |

496 | for i in maps: |

497 | points = [] |

498 | for j in i: |

499 | points += j |

500 | if g is None: |

501 | g = line(j,rgbcolor=c2) |

502 | else: |

503 | g += line(j,rgbcolor=c2) |

504 | if f is None: |

505 | f = polygon(points,rgbcolor=c1) |

506 | else: |

507 | f += polygon(points,rgbcolor=c1) |

508 | f += g |

509 | return f.show(axes=False) |

510 | |

511 | |

512 | |

513 | |

514 |