Search code examples
pythonnumpyscikit-imageangle

Python - is there a way to identify the directionality of angles?


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! The Teardrop Molecule. I want to measure the angles between subsequent segments of the trace in red.


Solution

  • 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
    

    enter image description here

    The sign change now indicates where the curve changes direction.