in the image below, a so-called teardrop-like DNA molecule is shown. For me, the distribution of the bending angles of this DNA is of very high interest. In red, a trace yielded by skimage.morphology.skeletonize
is shown.
In my program, the trace is being represented as a numpy-array: Each entry in the numpy-array is a float64 tuple corresponding to the x and y coordinate of the trace point in the image. So if I apply np.diff(trace, axis=0)
, I get the segments I want to analyze in vector from.
The problem: Now, I now how to calculate an angle between two subsequent segments. However, I do not now if the angle is oriented clockwise or counter-clockwise. For example, at one point, the DNA bends to the left, so the angle there would be different from the one measured when to DNA bends to the right, which it mostly does. What would be a way to fix this?
My idea so far: I thought of interpolating each segment and then checking if the next segment lies below or above the interpolated line, indicating a clockwise vs. counter-clockwise angle. But this seems hard to implement, and since I am a beginner, any input would be highly appreciated!
I think the cross product np.cross
will help you here. For a simple triangle going around first counter-clockwise, then clockwise this could look like this (can be optimized, left it longer for clarity). The angle output is the angle enclosed between the vector with the directions following the right-hand rule: EDIT: added a curve with direction change, added graph (annotations for going counter-clockwise):
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path
def drangl(vsgi, vsgi1):
'''sin(theta)=|axb|/(|a||b|)'''
return np.arcsin(np.cross(vsgi, vsgi1)/np.linalg.norm(vsgi)/np.linalg.norm(vsgi1))
print('going around the traingle counter-clockwise and clockwise:\n')
trace = np.array([(0.,0.), (10.,-1.), (5.,10.), (0.,0.)])
vsegs = np.diff(trace, axis=0)
for iv in [0,1,-1]: # counter-clockwise
print(f'ctr-cl: {iv:d} {drangl(vsegs[iv], vsegs[iv+1]):2.2f}')
for iv in [-1,1,0]: # clockwise
print(f'cw: {iv:d} {drangl(vsegs[iv], vsegs[iv+1]):2.2f}')
print('\nnow with direction change in the path:\n')
trace = np.array([(0.,0.), (10.,-1.), (9., 5.), (12., 5.), (5.,10.), (0.,0.)])
vsegs = np.diff(trace, axis=0)
for iv in [0,1,2,3,-1]: # counter-clockwise
print(f'ctr-cl: {iv:d} {drangl(vsegs[iv], vsegs[iv+1]):2.2f}')
for iv in [-1,3,2,1,0]: # clockwise
print(f'cw: {iv:d} {drangl(vsegs[iv], vsegs[iv+1]):2.2f}')
codes = [
Path.MOVETO,
Path.LINETO,
Path.LINETO,
Path.LINETO,
Path.LINETO,
Path.CLOSEPOLY,
]
path = Path(trace, codes)
fig, ax = plt.subplots()
patch = patches.PathPatch(path, facecolor='white', lw=2)
ax.add_patch(patch)
ax.set_xlim(-1, 13)
ax.set_ylim(-2, 11)
ax.annotate('turns right', xy=(9, 5), xytext=(5, 5),
arrowprops=dict(facecolor='black', shrink=0.05),
)
ax.annotate('turns left', xy=(12, 5), xytext=(10, 10),
arrowprops=dict(facecolor='black', shrink=0.05),
)
plt.show()
which produces
going around the traingle counter-clockwise and clockwise:
ctr-cl: 0 1.04
ctr-cl: 1 0.89
ctr-cl: -1 1.21
cw: -1 1.21
cw: 1 0.89
cw: 0 1.04
now with direction change in the path:
ctr-cl: 0 1.31
ctr-cl: 1 -1.41
ctr-cl: 2 0.62
ctr-cl: 3 1.41
ctr-cl: -1 1.21
cw: -1 1.21
cw: 3 1.41
cw: 2 0.62
cw: 1 -1.41
cw: 0 1.31
The sign change now indicates where the curve changes direction.