diff --git a/dnarbd/segmentmodel.py b/dnarbd/segmentmodel.py index cf82cd37adfd29d24f862e58abe1fe3627ba7104..cf834d67a93cff22d93eaf4daa14a082f4e13bfe 100644 --- a/dnarbd/segmentmodel.py +++ b/dnarbd/segmentmodel.py @@ -1410,8 +1410,9 @@ class SegmentModel(ArbdModel): filter_fn = lambda x: x is not None and x not in beads bs = list( filter( filter_fn, B.particle.intrahelical_neighbors ) ) beads.update(bs) - bs = list( filter( filter_fn, [n for b in bs for n in b.intrahelical_neighbors] ) ) - beads.update(bs) + for i in range(3): + bs = list( filter( filter_fn, [n for b in bs for n in b.intrahelical_neighbors] ) ) + beads.update(bs) beads = list(beads) @@ -1433,12 +1434,21 @@ class SegmentModel(ArbdModel): beads.remove(b) contours = [b.get_contour_position(s) for b in beads] + contours = np.array(contours, dtype=np.float16) # deliberately use low precision contours,ids = np.unique(contours, return_index=True) - if np.any( (contours[:-1] - contours[1:])**2 < 1e-6 ): + if np.any( (contours[:-1] - contours[1:])**2 < 1e-8 ): pdb.set_trace() ## TODO: keep closest beads beyond +-1.5 if there are fewer than 2 beads - tmp = [ci for ci in zip(contours,ids) if np.abs(ci[0]-0.5) < 1] + tmp = [] + dist = 1 + while len(tmp) < 5 and dist < 3: + tmp = [ci for ci in zip(contours,ids) if np.abs(ci[0]-0.5) < dist] + dist += 0.1 + + if len(tmp) <= 1: + raise Exception("Failed to fit spline into segment {}".format(s)) + contours = [ci[0] for ci in tmp] ids = [ci[1] for ci in tmp] @@ -1514,11 +1524,8 @@ class SegmentModel(ArbdModel): # pdb.set_trace() quats = np.array(quats) - ## TODOTODO test smoothing - try: - tck, u = interpolate.splprep( quats.T, u=contours, s=3, k=3 ) - except: - tck, u = interpolate.splprep( quats.T, u=contours, s=0, k=1 ) + # tck, u = interpolate.splprep( quats.T, u=contours, s=3, k=3 ) ;# cubic spline not as good + tck, u = interpolate.splprep( quats.T, u=contours, s=0, k=1 ) s.quaternion_spline_params = tck def _generate_bead_model(self,