两个函数的交点,多个x轴和y轴坐标

Intersection points between two functions, multiple x-axis and y-axis coordinates

我正在编写一个脚本,该脚本从 MS Excel 工作簿读取数据并进行一些绘图等操作。从 Excel 读取的数据是a_xa_ya_z 方向的加速度测量和 s 中的时间(单独的 numpy 数组)。在绘制之前,首先使用 5Hz 低通滤波器 过滤加速度,请参见图 1 的绘图示例,这是允许的加速度水平。

我需要找出超出限制曲线的时间量,但情节与 a_zabs(a_y) 有关,而不是时间。我尝试的解决方案是:

  1. 求加速度与极限曲线的交点
  2. 找到离交点最近的数据点(a_zabs(a_y))。
  3. 获取最接近交点的数据点的索引。
  4. 使用找到的时间数组索引并将两者相减以获得超出限制曲线的时间量。

问题是我的子弹 3 失败了。我已经成功地找到了限制曲线和过滤数据之间的交点。我也设法找到了最接近交点的数据点,但是,如果我找到了 a_z 最近的交点,这与我从 abs(a_y).[=36= 获得的索引不匹配]

通过以下方式找到交点:

f_l1 = LineString(np.column_stack((x, y1)))
s_l1 = LineString(np.column_stack((az3, abs(ay3))))
intersection1 = f_l1.intersection(s_l1)
az3_1,ay3_1 = LineString(intersection1).xy

所以我使用作为 LineString 导入的 shapely.geometry 来查找极限曲线(这里显示的是极限曲线 y1)和函数 s_l1(az,abs(a_y)).

为了找到最接近交点的数据点,我使用了以下方法:

我用来获取最接近交点的数据点的函数是:

def find_nearest_index(array, value):
    array = np.asarray(array)
    idx = np.abs(abs(array-value)).argmin()
    return idx

,其中数组是我过滤后的加速度数组,值是交点的 a_za_y 值。

我怀疑取决于 a_zabs(a_y) 的不同索引的原因是因为我的数据点与实际相交坐标“太远”,即我可能会得到 a_z 坐标,该值靠近交点但对应 abs(a_y) 太远了。 abs(a_y) 交集 point/data-point 相关性存在类似问题。对于即将进行的测量,我会增加采样频率,但我怀疑这能否完全解决问题。

我尝试了一些不同的方法但运气不佳,我的最新想法是在定位最近的数据点时尝试使用两个交点,所以我检查我从 find_nearest_index-function 获得的索引是否使用a_z 与我对 abs(a_y) 使用 find_nearest_index-function 得到的索引相同,但我不知道 how/if 是否可行。但也许我的问题有一个明显的解决方案,但我只是没有看到。

加速度的正确解决方案如下,其中我的数据点的索引与交点相匹配: Desirable plot between intersection points 然后,这些指数用于通过取 Delta_t=t[index2]-t[index1].

计算超出限制曲线的时间量

但是我通常会得到类似的东西,其中使用 a_z 找到的索引与使用 a_y) 找到的索引不同,导致错误的绘图,因此也是错误的 Delta_t : Typical plot between intersection points

这是我重用np.diff()的想法的方法。它允许轻松分段,因此可以获得所需的时间戳。小的修改将允许递归使用,因此在三个限制曲线的情况下应用简单。

import matplotlib.pyplot as plt
import numpy as np

### this is somewhat one of the bounds given in the OP
pathx = np.array([
    -1.7,
    -1.5, -0.5,
    0, 
    1.75, 5.4,
    6
])

pathy = np.array([
    0,
    0.75, 2.25,
    2.45,
    2.2, 0.75, 0
])

path = np.column_stack( ( pathx, pathy ) )

### this creates a random path
def rand_path( n ):
    vy = 0.04
    vx = 0
    xl = [0]
    yl = [0]
    for i in range( n ):
        phi = (1-1.6 *np.random.random() ) * 0.1
        mx = np.array(
            [
                [ +np.cos( phi ), np.sin( phi ) ],
                [ -np.sin( phi ), np.cos( phi ) ]
            ]
        )
        v = np.dot( mx, np.array( [ vx, vy ] ) )
        vx = v[0]
        vy = v[1]
        x = xl[-1] + vx
        y = yl[-1] + vy
        xl = xl + [ x ]
        yl = yl + [ y ]
    return xl, np.abs( yl )

### my version to check inside or not
def check_inside( point, path ):
    """
    check if point is inside convex boundary
    it is based on the sign of a cross product
    """
    out = 1
    for p2, p1 in zip( path[ 1: ], path[ :-1 ] ):
        q = p2 - p1
        Q = point - p1
        cross = q[0] * Q[1] - q[1] * Q[0]
        if cross > 0:
            out = 0
            break
    return out

### test data
xl ,yl = rand_path( 900 )
data = np.column_stack( ( xl, yl ) )

##check and use np.diff like in the other example 
cp = np.fromiter(
    ( check_inside( p, path ) for p in data ),
    int
)
ip = np.argwhere( np.diff( cp ) )

### get the points
if len( ip ):
    ip = np.concatenate( ip )
ipout = list()
for p in ip:
    if cp[ p ]:
        ipout.append( p + 1 )
    else:
        ipout.append( p )
        
pnts = data[ ipout ]

### split the line segments
innerSegments = list()
outerSegments = list()
ipchecklist= [0] + ipout + [ len( cp ) - 2 ]

for cntr,(s,e) in enumerate(
    zip( ipchecklist[:-1], ipchecklist[1:] )
):
    if cntr % 2:
        ss = s
        ee = e + 1
    else:
        ss = s + 1
        ee = e
    segment = data[ ss : ee ]
    if cntr % 2:
        outerSegments.append( segment )
    else:
        innerSegments.append( segment )

"""
Here would have only the points that are truly outside the border. 
Instead of getting the line segment data we could access a corresponding
time stamp in the same way and calculate the time outside the limit.
Having the points all outside would result in a systematic error towards
smaller times. We could also interpolate the crossing and the according
time to estimate the time of the crossing. This would remove the bias 
and improve the estimate.

With more than one boundary, we could feed the outside segments in the
same type of algorithm and repeat the procedure with a different border.
We all is put in a function, this would be a sort of recursive procedure.
"""
### plot
fig = plt.figure()

ax = fig.add_subplot( 1, 1, 1)
ax.scatter( pnts[:,0], pnts[:,1])

ax.plot( path[:,0], path[:,1])
# ~ ax.plot( xl, yl, ls='', marker='+')
for cnt, s in enumerate( innerSegments ):
    col= "#{:02x}0000".format(
        int( 25 + 230 * cnt / len( innerSegments ) )
    )
    ax.plot( s[:,0], s[:,1], color=col )
for cnt, s in enumerate( outerSegments ):
    col= "#0000{:02x}".format(
        int( 25 + 230 * (1 - cnt / len( outerSegments ) ) )
    )
    ax.plot( s[:,0], s[:,1], color=col )

ax.grid()
plt.show()