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函数是利用内部点的二阶精确中心差和边界上的一阶或二阶精确单面(正反向)差来计算的，我确信有一个函数可以逆转这一过程。

我在想，把原始点(t=0 at x1，y1)到Vx，Vy平面上的任意一点(xi，yi)之间的线导数，会给出速度分量之和。然后，我可以将这个值除以两点之间的距离，得到所需的时间。

这种方法有效吗？如果是这样的话，哪个numpy集成函数将是最好的应用？

​

​

我的数据的一个例子可以在这里找到http://www.filedropper.com/calculatearrivaltimefromgradientvalues060820

您的帮助将不胜感激。

编辑：

也许这张简化的图画也许能帮助我理解我想要达到的目的。

​

编辑：

感谢@一个与这段代码相连的人。我试着用间距为0.5x0.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')

可以在这里找到Inputdata，链接

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;DR；

在这个问题上，您需要应对多种挑战，主要是：

• 从梯度(矢量场)重建势场(标量场)

但也包括：

• 非矩形网格凹形船体的观测；
• 数值二维线积分和数值不精确；

它似乎可以通过选择一种特殊的插值和一种聪明的集成方式来解决(正如@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)

注意，这个内插器的接口稍有不同(注意如何传递parmaters )。

结果如下：

​

​

我们可以看到凸包外的区域可以外推(RBF是无网格的)。因此，选择adhoc插值绝对是解决问题的关键。但我们仍然需要意识到，外推法可能表现良好，但在某种程度上是毫无意义和危险的。

解决你的问题

@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]])

用经典的线性插值。结果是非常正确的，但我们在左下角仍然有一个边缘效应：

​

​

正如您注意到的，这个问题发生在集成开始并缺乏参考点的区域的轴中部。

Answer 2

以下是一种方法：

首先，为了能够进行集成，在常规网格上运行是很好的。在这里使用变量名x和y作为triang.x和triang.y的缩写，我们可以首先创建一个网格：

import numpy as np
n = 200 # Grid density
stepx = (max(x) - min(x)) / n
stepy = (max(y) - min(y)) / n
xspace = np.arange(min(x), max(x), stepx)
yspace = np.arange(min(y), max(y), stepy)
xgrid, ygrid = np.meshgrid(xspace, yspace)

然后，我们可以使用相同的dx函数在网格上插值dy和LinearTriInterpolator：

fdx = LinearTriInterpolator(masked_triang, dx)
fdy = LinearTriInterpolator(masked_triang, dy)

dxgrid = fdx(xgrid, ygrid)
dygrid = fdy(xgrid, ygrid)

现在来了整合的部分。原则上，我们选择的任何路径都应该使我们达到相同的价值。在实践中，由于存在缺失值和不同的密度，路径的选择对于得到一个合理准确的答案是非常重要的。

下面，我选择在n/2处沿x方向(从0到网格中部)在dxgrid上进行积分，然后在y方向上从0到i兴趣点在dygrid上进行积分。然后再从n/2到j点对dxgrid进行遍历。这是一种简单的方法，可以通过简单地选择主要位于数据范围“中间”的路径来确保大多数集成路径都在大量可用数据中。其他替代考虑将导致不同的路径选择。

所以我们会：

dxintegral = np.nancumsum(dxgrid, axis=1) * stepx
dyintegral = np.nancumsum(dygrid, axis=0) * stepy

然后(为了清晰起见，用某种蛮力)：

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

valintegral的结果是一个任意的常量，它可以帮助把“零”放在你想要的地方。

您的数据显示在这里：

ax.tricontourf(masked_triang, time_array)

​

这就是我在使用这种方法时要重建的内容：

ax.contourf(xgrid, ygrid, valintegral)

​

希望这会有所帮助。

如果要重新访问原始三角剖分点上的值，可以在valintegral常规网格数据上使用valintegral。

编辑：

在回复您的编辑时，您的上述修改有以下几个错误：

1. 将行(dx,dy) = np.gradient(grid_z1)更改为(dy,dx) = np.gradient(grid_z1)
2. 在集成循环中，将dyintegral[i, len(yy) // 2]项更改为dyintegral[i, len(xx) // 2]
3. 最好将行valintegral = valintegral * np.isfinite(dxintegral)替换为valintegral[np.isnan(dx)] = np.nan