-
Notifications
You must be signed in to change notification settings - Fork 8
/
dna
executable file
·162 lines (142 loc) · 4.69 KB
/
dna
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
#!/usr/bin/env python
from __future__ import division,absolute_import
from geode import *
from geode.geometry.platonic import *
from gui import *
import subprocess
import sys
# Useful commands:
#
# convert -crop 211x112+740+231 +repage -delay 4 -loop 0 `seq -f 'capture-%03g.png' 0 149` dna-large.gif
# convert -crop 211x112+740+231 +repage -resize 32x17 -delay 4 -loop 0 `seq -f 'capture-%03g.png' 0 149` dna-small.gif
# Properties
props = PropManager()
z_scale = props.add('z_scale',.5).set_category('dna')
backbone_width = props.add('backbone_width',.15).set_category('dna')
base_pair_width = props.add('base_pair_width',.1).set_category('dna')
periods = props.add('periods',2).set_category('dna')
resolution = props.add('resolution',50).set_category('dna')
spacing = props.add('spacing',10).set_category('dna')
around = props.add('around',10).set_category('dna')
capture_frames = props.add('capture_frames',150+10).set_category('dna')
seed = props.add('seed',31833).set_category('dna')
mode = props.add('mode','cylinder').set_allowed('flat cylinder'.split()).set_category('dna')
# Geometry of B-DNA
# For details, see http://iopscience.iop.org/1367-2630/15/9/093008/pdf/1367-2630_15_9_093008.pdf
pitch = pi/180*38
H_bar = tan(pitch) # tan pitch = H/(2*pi*a), H_bar = H/(2*pi)
phase = pi/180*140
bases_per_turn = 10.5
bases = 'atgc'
# http://wiki.answers.com/Q/What_is_the_color_of_the_base_pairs_in_DNA
if 0: # Jellybean colors
base_colors = ((0,0,1),(1,1,0),(0,1,0),(1,0,0))
backbone_colors = ((.75,0,.25),(.25,0,.75))
else:
base_colors = ((.2,.1,.9),(.9,.8,0),(0,.8,.2),(.9,.2,.1))
backbone_colors = ((.9,.1,.25),(.4,.1,.55))
backbone_colors = ((.75,0,.25),(.25,0,.75))
@cache
def sequence():
n = 100
random.seed(seed())
return random.randint(4,size=n)
@cache
def helix_curve():
n = periods()
t = linspace(0,2*pi*n,num=resolution()*n)
X = empty((len(t),3))
s = time()
X[:,0] = cos(t+s)
X[:,1] = sin(t+s)
X[:,2] = H_bar*t
return X
def span_start(i):
s = time()
t = 2*pi/bases_per_turn*i
return asarray([cos(t+s),sin(t+s),H_bar*t])
@cache
def helix_mesh():
core = helix_curve()
if mode()=='cylinder':
return cylinder_topology(len(core)-1,around())
elif mode()=='flat':
return grid_topology(len(core)-1,1)
@cache
def helix_X():
core = helix_curve()
if mode()=='cylinder':
return revolve_around_curve(core,backbone_width(),around())[1]
elif mode()=='flat':
tangent = core[1:]-core[:-1]
tangent = concatenate([[tangent[0]],(tangent[:-1]+tangent[1:])/2,[tangent[-1]]])
assert tangent.shape==core.shape
outwards = core.copy()
outwards[:,2] = 0
outwards = normalized(outwards)
left = normalized(cross(outwards,tangent))
X = empty((len(core),2,3))
shift = backbone_width()*left
X[:,0] = core-shift
X[:,1] = core+shift
return X.reshape(-1,3)
@cache
def other_rotation():
return Rotation.from_angle_axis(phase,(0,0,1))
@cache
def other_helix_X():
return other_rotation()*helix_X()
def span(i,j):
def span():
x0 = span_start(i)
x0[:2] *= .99
x1 = other_rotation()*x0
mid = (x0+x1)/2
R = base_pair_width()
return open_cylinder_mesh(mid,(x0,x1)[j],R,around())
return cache(span)
def span_scene(i,j):
sp = span(i,j)
color = base_colors[sequence()[i]^j]
return MeshScene(props,cache(lambda:sp()[0]),cache(lambda:sp()[1]),color,color)
# View
app = QEApp(sys.argv,True)
main = MainWindow(props)
main.resize_timeline(80)
main.timeline.info.set_loop(1)
main.view.clear_color = zeros(3)
props.get('last_frame').set(1000000)
frame = props.get('frame')
frame_rate = props.get('frame_rate')
# Define time so that 150 frames is 2pi
time = cache(lambda:2*pi/150*frame())
# Key bindings
main.add_menu_item('Timeline','Play/Stop',main.timeline.info.set_play,'Ctrl+p')
main.add_menu_item('Timeline','Step back',main.timeline.info.go_back,'Ctrl+< Ctrl+,')
main.add_menu_item('Timeline','Step forward',main.timeline.info.go_forward,'Ctrl+.') # Possibly due to a Qt bug, 'Ctrl+> Ctrl+.' doesn't work
# Video capture
def capture():
print 'capturing'
n = capture_frames()
i = [0]
def grab():
filename = 'capture-%03d.png'%i[0]
print ' writing %s'%filename
subprocess.check_output(['scrot','-u',filename])
i[0] += 1
if i[0]==n:
main.view.post_render = None
main.timeline.info.set_play()
main.view.post_render = grab
main.timeline.info.set_play()
main.add_menu_item('Capture','Capture',capture,'Ctrl+g')
# Add scenes
for j in 0,1:
color = backbone_colors[j]
main.view.add_scene('helix %d'%j,MeshScene(props,helix_mesh,(helix_X,other_helix_X)[j],color,color))
for i in xrange(int(bases_per_turn*periods())):
for j in 0,1:
main.view.add_scene('span %d %d'%(i,j),span_scene(i,j))
# Launch!
main.init()
app.run()