numpy.gradient 函数的反函数

Question

我需要创建一个与 np.gradient 函数相反的函数。

其中 Vx,Vy 数组（速度分量向量）是输入，输出将是数据点 x,y 处的反导数（到达时间）数组。

我在 (x,y) 网格上有数据，每个点都有标量值（时间）。

我使用了 numpy 梯度函数和线性插值来确定每个点的梯度向量速度（Vx，Vy）（见下文）。

我通过以下方式实现了这一目标：

 #LinearTriInterpolator applied to a delaunay triangular mesh
 LTI= LinearTriInterpolator(masked_triang, time_array)

 #Gradient requested at the mesh nodes:
 (Vx, Vy) = LTI.gradient(triang.x, triang.y)

下面的第一幅图显示了每个点的速度向量，点标签表示形成导数 (Vx,Vy) 的时间值

下图显示了导数 (Vx,Vy)的结果标量值，绘制为带有关联节点标签的彩色等高线图。

所以我的挑战是：

我需要逆转这个过程！

使用梯度向量 (Vx,Vy) 或结果标量值来确定该点的原始时间值。

这可能吗？

知道 numpy.gradient 函数是使用内部点的二阶精确中心差异和边界处的一阶或二阶精确一侧（向前或向后）差异计算的，我相信有一个函数可以逆转这一点过程。

我在想，在原始点（x1,y1 处的 t=0）到 Vx,Vy 平面上的任何点 (xi,yi) 之间取线导数会给我速度分量的总和。然后我可以用这个值除以两点之间的距离来得到所用的时间..

这种方法行得通吗？如果是这样，最好应用哪个 numpy 集成函数？

可以在此处找到我的数据示例 [http://www.filedropper.com/calculatearrivaltimefromgradientvalues060820]

您的帮助将不胜感激

编辑：

也许这个简化的绘图可能有助于理解我想要到达的地方..

编辑：

感谢@Aguy 提供了此代码。我尝试使用间距为 0.5 x 0.5m 的网格并计算每个网格点的梯度来获得更准确的表示，但是我无法正确集成它。我也有一些影响我不知道如何纠正的结果的边缘影响。

import numpy as np
from scipy import interpolate
from matplotlib import pyplot
from mpl_toolkits.mplot3d import Axes3D

#Createmesh grid with a spacing of 0.5 x 0.5
stepx = 0.5
stepy = 0.5
xx = np.arange(min(x), max(x), stepx)
yy = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xx, yy)
grid_z1 = interpolate.griddata((x,y), Arrival_Time, (xgrid, ygrid), method='linear') #Interpolating the Time values

#Formatdata
X = np.ravel(xgrid)
Y= np.ravel(ygrid)
zs = np.ravel(grid_z1)
Z = zs.reshape(X.shape)

#Calculate Gradient
(dx,dy) = np.gradient(grid_z1) #Find gradient for points on meshgrid

Velocity_dx= dx/stepx #velocity ms/m
Velocity_dy= dy/stepx #velocity ms/m

Resultant = (Velocity_dx**2 + Velocity_dy**2)**0.5 #Resultant scalar value ms/m

Resultant = np.ravel(Resultant)

#Plot Original Data F(X,Y) on the meshgrid
fig = pyplot.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x,y,Arrival_Time,color='r')
ax.plot_trisurf(X, Y, Z)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.set_zlabel('Time (ms)')
pyplot.show()

#Plot the Derivative of f'(X,Y) on the meshgrid
fig = pyplot.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(X,Y,Resultant,color='r',s=0.2)
ax.plot_trisurf(X, Y, Resultant)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.set_zlabel('Velocity (ms/m)')
pyplot.show()

#Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1)*stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0)*stepy

valintegral = np.ma.zeros(dxintegral.shape)
for i in range(len(yy)):
    for j in range(len(xx)):
        valintegral[i, j] = np.ma.sum([dxintegral[0, len(xx) // 2], 
    dyintegral[i, len(yy)  // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral = valintegral * np.isfinite(dxintegral)

现在 np.gradient 应用于每个网格节点 (dx,dy) = np.gradient(grid_z1)

现在在我的过程中，我将分析上面的梯度值并进行一些调整（正在创建一些我需要纠正的不寻常的边缘效果）然后将这些值整合到一个非常类似于f(x,y) 如上所示。

我需要一些帮助来调整积分功能：

#Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1)*stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0)*stepy

valintegral = np.ma.zeros(dxintegral.shape)
for i in range(len(yy)):
    for j in range(len(xx)):
        valintegral[i, j] = np.ma.sum([dxintegral[0, len(xx) // 2], 
    dyintegral[i, len(yy)  // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral = valintegral * np.isfinite(dxintegral)

现在我需要在原始 (x,y) 点位置计算新的“时间”值。

更新（08-09-20）：在@Aguy 的帮助下，我得到了一些有希望的结果。结果如下所示（蓝色轮廓代表原始数据，红色轮廓代表积分值）。

我仍在研究一种集成方法，该方法可以消除 min(y) 和 max(y) 区域的不准确性

from matplotlib.tri import (Triangulation, UniformTriRefiner, 
CubicTriInterpolator,LinearTriInterpolator,TriInterpolator,TriAnalyzer)
import pandas as pd
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

#-------------------------------------------------------------------------
# STEP 1: Import data from Excel file, and set variables
#-------------------------------------------------------------------------
df_initial = pd.read_excel(
r'C:\Users\morga\PycharmProjects\venv\Development\Trial'
r'.xlsx')

输入数据可以在这里找到链接

df_initial = df_initial .sort_values(by='Delay', ascending=True) #Update dataframe and sort by Delay
x = df_initial ['X'].to_numpy() 
y = df_initial ['Y'].to_numpy() 
Arrival_Time = df_initial ['Delay'].to_numpy() 

# Createmesh grid with a spacing of 0.5 x 0.5
stepx = 0.5
stepy = 0.5
xx = np.arange(min(x), max(x), stepx)
yy = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xx, yy)
grid_z1 = interpolate.griddata((x, y), Arrival_Time, (xgrid, ygrid), method='linear')  # Interpolating the Time values

# Calculate Gradient (velocity ms/m)
(dy, dx) = np.gradient(grid_z1)  # Find gradient for points on meshgrid


Velocity_dx = dx / stepx  # x velocity component ms/m
Velocity_dy = dy / stepx  # y velocity component ms/m

# Integrate to compare the original data input
dxintegral = np.nancumsum(Velocity_dx, axis=1) * stepx
dyintegral = np.nancumsum(Velocity_dy, axis=0) * stepy

valintegral = np.ma.zeros(dxintegral.shape)  # Makes an array filled with 0's the same shape as dx integral
for i in range(len(yy)):
    for j in range(len(xx)):
        valintegral[i, j] = np.ma.sum(
        [dxintegral[0, len(xx) // 2], dyintegral[i, len(xx) // 2], dxintegral[i, j], - dxintegral[i, len(xx) // 2]])
valintegral[np.isnan(dx)] = np.nan
min_value = np.nanmin(valintegral)

valintegral = valintegral + (min_value * -1)

##Plot Results

fig = plt.figure()
ax = fig.add_subplot()
ax.scatter(x, y, color='black', s=7, zorder=3)
ax.set_xlabel('X-Coordinates')
ax.set_ylabel('Y-Coordinates')
ax.contour(xgrid, ygrid, valintegral, levels=50, colors='red', zorder=2)
ax.contour(xgrid, ygrid, grid_z1, levels=50, colors='blue', zorder=1)
ax.set_aspect('equal')
plt.show()

Answer 1

TL; 博士;

在这个问题上你有很多挑战需要解决，主要是：

• 从其梯度（矢量场）进行电位重建（标量场）

但是也：

• 在具有非矩形网格的凹壳中观察；
• 数值二维线积分和数值不准确；

似乎可以通过选择临时插值和智能集成方式来解决（如 所指出的@Aguy）。

MCVE

第一次，让我们构建一个 MCVE 来突出上述关键点。

数据集

我们重新创建一个标量场及其梯度。

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt

def f(x, y):
    return x**2 + x*y + 2*y + 1

Nx, Ny = 21, 17
xl = np.linspace(-3, 3, Nx)
yl = np.linspace(-2, 2, Ny)

X, Y = np.meshgrid(xl, yl)
Z = f(X, Y)
zl = np.arange(np.floor(Z.min()), np.ceil(Z.max())+1, 2)

dZdy, dZdx = np.gradient(Z, yl, xl, edge_order=1)
V = np.hypot(dZdx, dZdy)

标量字段如下所示：

axe = plt.axes(projection='3d')
axe.plot_surface(X, Y, Z, cmap='jet', alpha=0.5)
axe.view_init(elev=25, azim=-45)

而且，向量场看起来像：

axe = plt.contour(X, Y, Z, zl, cmap='jet')
axe.axes.quiver(X, Y, dZdx, dZdy, V, units='x', pivot='tip', cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()

事实上，梯度对于潜在水平是正常的。我们还绘制了梯度幅度：

axe = plt.contour(X, Y, V, 10, cmap='jet')
axe.axes.set_aspect('equal')
axe.axes.grid()

原始场重建

如果我们天真地从梯度重建标量场：

SdZx = np.cumsum(dZdx, axis=1)*np.diff(xl)[0]
SdZy = np.cumsum(dZdy, axis=0)*np.diff(yl)[0]

Zhat = np.zeros(SdZx.shape)
for i in range(Zhat.shape[0]):
    for j in range(Zhat.shape[1]):
        Zhat[i,j] += np.sum([SdZy[i,0], -SdZy[0,0], SdZx[i,j], -SdZx[i,0]])
        
Zhat += Z[0,0] - Zhat[0,0]

我们可以看到全局结果大致正确，但在梯度幅度较低的情况下，级别不太准确：

插值场重建

如果我们增加网格分辨率并选择一个特定的插值（通常在处理网格时），我们可以获得更精细的场重建：

r = np.stack([X.ravel(), Y.ravel()]).T
Sx = interpolate.CloughTocher2DInterpolator(r, dZdx.ravel())
Sy = interpolate.CloughTocher2DInterpolator(r, dZdy.ravel())

Nx, Ny = 200, 200
xli = np.linspace(xl.min(), xl.max(), Nx)
yli = np.linspace(yl.min(), yl.max(), Nx)
Xi, Yi = np.meshgrid(xli, yli)
ri = np.stack([Xi.ravel(), Yi.ravel()]).T

dZdxi = Sx(ri).reshape(Xi.shape)
dZdyi = Sy(ri).reshape(Xi.shape)

SdZxi = np.cumsum(dZdxi, axis=1)*np.diff(xli)[0]
SdZyi = np.cumsum(dZdyi, axis=0)*np.diff(yli)[0]

Zhati = np.zeros(SdZxi.shape)
for i in range(Zhati.shape[0]):
    for j in range(Zhati.shape[1]):
        Zhati[i,j] += np.sum([SdZyi[i,0], -SdZyi[0,0], SdZxi[i,j], -SdZxi[i,0]])
        
Zhati += Z[0,0] - Zhati[0,0]

哪个绝对表现得更好：

所以基本上，使用临时插值器增加网格分辨率可能会帮助您获得更准确的结果。插值器还解决了从三角形网格中获取规则矩形网格以执行积分的需要。

凹凸包

您还指出了边缘的不准确性。这些是插值选择和集成方法相结合的结果。当标量场到达具有很少插值点的凹面区域时，积分方法无法正确计算标量场。在选择无电气嵌入时，问题消失了。

为了说明这一点，让我们从 MCVE 中删除一些数据：

q = np.full(dZdx.shape, False)
q[0:6,5:11] = True
q[-6:,-6:] = True
dZdx[q] = np.nan
dZdy[q] = np.nan

然后可以按如下方式构造插值：

q2 = ~np.isnan(dZdx.ravel())
r = np.stack([X.ravel(), Y.ravel()]).T[q2,:]
Sx = interpolate.CloughTocher2DInterpolator(r, dZdx.ravel()[q2])
Sy = interpolate.CloughTocher2DInterpolator(r, dZdy.ravel()[q2])

执行积分我们看到，除了经典的边缘效应，我们在凹面区域（船体凹面的摆动点划线）中的准确值也较低，并且我们在凸面体之外没有数据，因为 Clough Tocher 是基于网格的插值：

Vl = np.arange(0, 11, 1)
axe = plt.contour(X, Y, np.hypot(dZdx, dZdy), Vl, cmap='jet')
axe.axes.contour(Xi, Yi, np.hypot(dZdxi, dZdyi), Vl, cmap='jet', linestyles='-.')
axe.axes.set_aspect('equal')
axe.axes.grid()

所以基本上我们在拐角处看到的错误很可能是由于积分问题加上仅限于凸包的插值。

为了克服这个问题，我们可以选择一个不同的插值，例如 RBF（径向基函数内核），它能够在凸包外创建数据：

Sx = interpolate.Rbf(r[:,0], r[:,1], dZdx.ravel()[q2], function='thin_plate')
Sy = interpolate.Rbf(r[:,0], r[:,1], dZdy.ravel()[q2], function='thin_plate')

dZdxi = Sx(ri[:,0], ri[:,1]).reshape(Xi.shape)
dZdyi = Sy(ri[:,0], ri[:,1]).reshape(Xi.shape)

注意这个插值器的接口略有不同（注意参数是如何传递的）。

结果如下：

我们可以看到可以外推凸包外的区域（RBF 是无网格的）。所以选择adhoc interpolant绝对是解决你问题的关键。但是我们仍然需要意识到外推可能表现良好，但在某种程度上毫无意义且危险。

解决您的问题

提供的答案@Aguy非常好，因为它设置了一种巧妙的积分方式，不受凸包外缺失点的干扰。但是正如您所提到的，凸包内的凹区域不准确。

如果您希望消除检测到的边缘效应，您将不得不求助于能够进行外推的插值器，或者寻找另一种集成方法。

插值变化

使用 RBF 插值似乎可以解决您的问题。这是完整的代码：

df = pd.read_excel('./Trial-Wireup 2.xlsx')
x = df['X'].to_numpy()
y = df['Y'].to_numpy()
z = df['Delay'].to_numpy()

r = np.stack([x, y]).T

#S = interpolate.CloughTocher2DInterpolator(r, z)
#S = interpolate.LinearNDInterpolator(r, z)
S = interpolate.Rbf(x, y, z, epsilon=0.1, function='thin_plate')

N = 200
xl = np.linspace(x.min(), x.max(), N)
yl = np.linspace(y.min(), y.max(), N)
X, Y = np.meshgrid(xl, yl)

#Zp = S(np.stack([X.ravel(), Y.ravel()]).T)
Zp = S(X.ravel(), Y.ravel())
Z = Zp.reshape(X.shape)

dZdy, dZdx = np.gradient(Z, yl, xl, edge_order=1)

SdZx = np.nancumsum(dZdx, axis=1)*np.diff(xl)[0]
SdZy = np.nancumsum(dZdy, axis=0)*np.diff(yl)[0]

Zhat = np.zeros(SdZx.shape)
for i in range(Zhat.shape[0]):
    for j in range(Zhat.shape[1]):
        #Zhat[i,j] += np.nansum([SdZy[i,0], -SdZy[0,0], SdZx[i,j], -SdZx[i,0]])
        Zhat[i,j] += np.nansum([SdZx[0,N//2], SdZy[i,N//2], SdZx[i,j], -SdZx[i,N//2]])
        
Zhat += Z[100,100] - Zhat[100,100]

lz = np.linspace(0, 5000, 20)
axe = plt.contour(X, Y, Z, lz, cmap='jet')
axe = plt.contour(X, Y, Zhat, lz, cmap='jet', linestyles=':')
axe.axes.plot(x, y, '.', markersize=1)
axe.axes.set_aspect('equal')
axe.axes.grid()

其图形呈现如下：

边缘效应消失了，因为 RBF 插值可以在整个网格上外推。您可以通过比较基于网格的插值的结果来确认。

线性

克拉夫·托彻

积分变量顺序变化

我们还可以尝试找到更好的方法来整合和减轻边缘效应，例如。让我们改变积分变量顺序：

Zhat[i,j] += np.nansum([SdZy[N//2,0], SdZx[N//2,j], SdZy[i,j], -SdZy[N//2,j]])

使用经典的线性插值。结果相当正确，但我们仍然在左下角有边缘效果：

正如您所注意到的，问题发生在积分开始且缺少参考点的区域中的轴的中间。