def METHOD_NAME(skeleton):
    axes = np.empty_like(skeleton, dtype=np.float32)
    axes[-1] = skeleton[-1] - skeleton[-2]
    axes[:-1] = skeleton[1:] - skeleton[:-1]
    axes[1:-1] = (axes[1:-1] + axes[:-2]) / 2 # use average of neighboring edge directions where available
    norm = np.linalg.norm(axes, axis=-1)
    axes[norm > 0] /= norm[norm > 0, None]
    return axes