Commit 3a7e63fc authored by cmaffeo2's avatar cmaffeo2
Browse files

Fix spline interpolation; added routines for transforming entire model

parent e6d4a78d
......@@ -150,7 +150,7 @@ def quaternion_from_matrix( R ):
def rotationAboutAxis(axis,angle, normalizeAxis=True):
if normalizeAxis: axis = axis / np.linalg.norm(axis)
angle *= 0.5 * np.pi/180
angle = angle * 0.5 * np.pi/180
cos = np.cos( angle )
sin = np.sin( angle )
q = [cos] + [sin*x for x in axis]
......
......@@ -420,7 +420,7 @@ class Segment(ConnectableElement, Group):
ids = list(filter(lambda i: position_filter(r[i,:]), ids))
return ids
def translate(self, translation_vector, about=None, position_filter=None, contour_filter=None):
def translate(self, translation_vector, position_filter=None, contour_filter=None):
dr = np.array(translation_vector)
tck, u = self.position_spline_params
r = self.contour_to_position(u)
......@@ -429,7 +429,7 @@ class Segment(ConnectableElement, Group):
## Translate
r[ids,:] = r[ids,:] + dr[np.newaxis,:]
self.set_position_splines(r,u)
self.set_splines(u,r)
def rotate(self, rotation_matrix, about=None, position_filter=None, contour_filter=None):
tck, u = self.position_spline_params
......@@ -445,16 +445,16 @@ class Segment(ConnectableElement, Group):
## TODO: do this more efficiently
r[ids,:] = np.array([rotation_matrix.dot(r[i,:]-dr) + dr for i in range(r.shape[0])])
self.set_position_splines(r,u)
self.set_splines(u,r)
if self.quaternion_spline_params is not None:
## TODO: performance: don't shift between quaternion and matrix representations so much
tck, u = self.quaternion_spline_params
orientations = [self.contour_to_orientation(v)) for v in u]
orientations = [self.contour_to_orientation(v) for v in u]
for i in ids:
orientations[i,:] = rotation_matrix.dot(orientations[i])
quats = [quaternion_from_matrix(o) for o in orientations]
self.set_orientation_splines(quats,u)
self.set_orientation_splines(u, quats)
def _set_splines_from_ends(self):
self.quaternion_spline_params = None
......@@ -1386,6 +1386,7 @@ class SegmentModel(ArbdModel):
self._bonded_potential = dict() # cache for bonded potentials
self._generate_strands()
self.potentials = []
self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist, escapable_twist)
self.useNonbondedScheme( nbDnaScheme )
......@@ -1540,6 +1541,26 @@ class SegmentModel(ArbdModel):
print("WARNING: called deprecated command '_update_segment_positions; use 'update_splines' instead")
return self.update_splines(bead_coordinates)
## Operations on spline coordinates
def translate(self, translation_vector, position_filter=None):
for s in self.segments:
s.translate(translation_vector, position_filter=position_filter)
def rotate(self, rotation_matrix, about=None, position_filter=None):
for s in self.segments:
s.rotate(rotation_matrix, about=about, position_filter=position_filter)
def get_center(self, include_ssdna=False):
if include_ssdna:
segments = self.segments
else:
segments = list(filter(lambda s: isinstance(s,DoubleStrandedSegment),
self.segments))
centers = [s.get_center() for s in segments]
weights = [s.num_nt*2 if isinstance(s,DoubleStrandedSegment) else s.num_nt for s in segments]
# centers,weights = [np.array(a) for a in (centers,weights)]
return np.average( centers, axis=0, weights=weights)
def update_splines(self, bead_coordinates):
""" Set new function for each segments functions
contour_to_position and contour_to_orientation """
......@@ -1626,14 +1647,8 @@ class SegmentModel(ArbdModel):
# print(" {}[{:03f}]: {:0.3f}".format(b.parent.name, b.contour_position, c) )
ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
tck = ret[0][0]
tck, u = interpolate.splprep( positions, u=contours, s=0, k=1 )
ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
tck2 = ret[0][0]
assert( tck[0][0] == tck2[0][0] )
# if len(beads) < 8:
# ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
# tck = ret[0][0]
......@@ -1652,7 +1667,7 @@ class SegmentModel(ArbdModel):
# if ret[2] > 0:
# pdb.set_trace()
s.position_spline_params = (tck,ret[1][0])
s.position_spline_params = (tck,u)
""" Get twist """
cb = [e for e in zip(contours,beads) if 'orientation_bead' in e[1].__dict__]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment