Ticket #11564: polyhedron_unfold.py

File polyhedron_unfold.py, 13.8 KB (added by QuantumKing, 8 years ago)
Line 
1from sage.geometry.polyhedra import cyclic_sort_vertices_2d
2from sage.plot.line import line
3from sage.plot.polygon import polygon
4from math import sqrt
5from copy import copy
6from random import random
7
8##########################################################
9#
10#     Returns the dot product between two vectors
11#
12##########################################################
13
14def 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
26def 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
38def 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   
50def 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
62def 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
72def 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
84def 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
149def 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
211def 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
245def 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
284def 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   
300def 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
363def 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
443def 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
473def 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