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

Question

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

我需要找到超出限制曲线的时间量，但该图与a_z 和abs(a_y) 相关，而不是时间。我尝试的解决方案是：

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

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

交点通过以下方式找到：

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)) 之间的交点。

为了找到最接近交叉点的数据点，我使用了以下方法：

Intersection of two graphs in Python, find the x value

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

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

，其中数组是我过滤后的加速度数组，值是交点的a_z 或a_y 值。

我怀疑不同索引取决于a_z 或abs(a_y) 的原因是因为我的数据点与实际交叉点坐标“太远”，即我可能会得到一个a_z 坐标，该值接近交点但对应的 abs(a_y) 很遥远。 abs(a_y) 交叉点点/数据点相关性也存在类似问题。对于即将进行的测量，我会提高采样频率，但我怀疑这会完全解决问题。

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

加速度的正确解决方案如下所示，其中我的数据点的索引与交点匹配： Desirable plot between intersection points 然后通过取 Delta_t=t[index2]-t[index1]，这些指数用于计算超出限制曲线的时间量。

但我通常会得到类似的结果，使用 a_z 找到的索引与使用 a_y) 找到的索引不同，导致绘图错误，因此 Delta_t 也错误： Typical plot between intersection points

Answer 1

这是我重用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()