gpt4 book ai didi

python - 基于网格的多元三次插值

转载 作者:太空宇宙 更新时间:2023-11-03 18:18:53 24 4
gpt4 key购买 nike

说我有兴趣

y = f(n, m, o)


我为 nmo创建离散化网格:

import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')


为了简单起见,我们假设 f(n,m,o) = n+m+p

Y = N+M+O


因此,我现在在 f内具有 Y的网格明智近似值。从我的角度来看 documentation


我需要看一下多元插值
scipy.interpolate.Griddata看起来不错,但它确实适用于非结构化网格。这意味着我必须将所有结果转换为长数组-这并没有利用我实际上拥有结构化网格的优势
结构化下列出的其他方法不能同时支持三次和多元(大于2维)输入。


给定我已经创建的变量,在此问题上进行插值的好方法是什么?

更新资料

这是第一个答案的实现:

import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')
Y = N+M+O

n,m,o = 5, 103, 1007

m_i = np.argmin(np.abs(np.floor(m)-mGrid))
m_f = m - m_i

n_i = np.argmin(np.abs(np.floor(n)-nGrid))
n_f = n - n_i

o_i = np.argmin(np.abs(np.floor(o)-oGrid))
o_f = o - o_i

A = Y[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]

# cubic smoothing polynome in 1D
# Catmull-Rom style
def v(p0, p1, p2, p3, f):
return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
p1 * (3*f**3 - 5*f**2 + 2) +
p2 * (-3*f**3 + 4*f**2 + f) +
p3 * (f**3 - f**2))



B = v(A[0], A[1], A[2], A[3], m_f)
C = v(B[0], B[1], B[2], B[3], n_f)
D = v(C[0], C[1], C[2], C[3], o_f)

# D is the final interpolated array

print m+n+o, D


不幸的是, f(n,m,o) = 1115,而 D = 2215。由于与该方法没有任何联系,因此我很难理解到底发生了什么以及为什么如此近似的原因。

最佳答案

如果找不到好东西,那就自己滚吧!

但请注意,三次三次插值有些棘手。下面描述简单的三次三次插值。根据您的需要,有更多更复杂,更快速的方法(基于矩阵)可用。还应该意识到以下事实:有几种三次外推方法,我只使用可能是最常见的一种(Catmull-Rom)。

在3D阵列M中查找点(m,n,o)时,其三次三次插值的周围由4x4x4点(或3x3x3单元)网格组成。这些点在每个轴上的定义是:每个轴的floor(a)-1floor(a)floor(a)+1floor(a)+2等。 (想想一个魔方,其中(m,n,o)点位于不可见的中心立方体中。)

然后三次三次插值结果是这64个点的加权平均值。权重取决于点到网格点的距离。

让我们定义:

m_i = np.floor(m).astype('int')    # integer part
m_f = m - m_i # fraction


以及对n和o的类似定义。现在我们需要对数组进行操作:

A = M[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]


现在,我们的点是(1 + m_f,1 + n_f,1 + o_f)。 (请注意,这假设矩阵是无限的。否则,必须考虑边缘效应。)

数组的64个系数可以用一些聪明的矩阵数学来计算,但是在这里足以知道插值是关联的,即我们一次可以做一个轴:

# cubic smoothing polynome in 1D
# Catmull-Rom style
def v(p0, p1, p2, p3, f):
return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
p1 * (3*f**3 - 5*f**2 + 2) +
p2 * (-3*f**3 + 4*f**2 + f) +
p3 * (f**3 - f**2))


# interpolate axis-by-axis
B = v(A[0], A[1], A[2], A[3], m_f)
C = v(B[0], B[1], B[2], B[3], n_f)
D = v(C[0], C[1], C[2], C[3], o_f)

# D is the final interpolated value


因此,插补一次只计算一个轴。



让我们将其变为更实际的东西。首先,我们需要从自己的坐标空间到数组坐标进行一些地址计算。如果定义了m轴:

m0: position of the first grid plane on the axis
m1: position of the last grid plane on the axis
n_m: number of samples


我们可以从“真实世界”位置 m计算等效数组位置:

m_arr = (n_m - 1) * (m -  m0) / (m1 - m0)


-1来自击剑现象(如果我们有10个样本,则第一个与最后一个样本之间的距离为9)。

现在我们可以将所有内容放到一个函数中:

# cubic smoothing polynome in 1D
# Catmull-Rom style
def interp_catmull_rom(p0, p1, p2, p3, f):
return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
p1 * (3*f**3 - 5*f**2 + 2) +
p2 * (-3*f**3 + 4*f**2 + f) +
p3 * (f**3 - f**2))

# linear interpolation
# only needs p1 and p2, the rest are there for compatibility
def interp_linear(p0, p1, p2, p3, f):
return (1-f)*p1 + f*p2


# don't interpolate, use the nearest point
# only needs p1 and p2, the rest are there for compatibility
def interp_nearest(p0, p1, p2, p3, f):
if f > .5:
return p2
else:
return p1

# 3D interpolation done axis-by-axis
def tri_interp(M, f, m, n, o, m0, m1, n0, n1, o0, o1):
# M: 3D array
# f: interpolation function to use
# m,n,o: coordinates where to interpolate
# m0: real world minimum for m
# m1: real world maximum for m
# n0,n1,o0,o0: as with m0 and m1

# calculate the array coordinates
m_arr = (M.shape[0] - 1) * (m - m0) / (m1 - m0)
n_arr = (M.shape[1] - 1) * (n - n0) / (n1 - n0)
o_arr = (M.shape[2] - 1) * (o - o0) / (o1 - o0)

# if we are out of our data box, return a nan
if m_arr < 0 or m_arr > M.shape[0] or \
n_arr < 0 or n_arr > M.shape[1] or \
o_arr < 0 or o_arr > M.shape[2]:
return np.nan

# calculate the integer parts and fractions
m_i = np.floor(m_arr).astype('int')
n_i = np.floor(n_arr).astype('int')
o_i = np.floor(o_arr).astype('int')
m_f = m_arr - m_i
n_f = n_arr - n_i
o_f = o_arr - o_i

# edge effects may be nasty, we may need elements outside of the array
# there may be more efficient ways to avoid it, but we'll create lists of
# coordinates:

n_coords, m_coords, o_coords = np.mgrid[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]

# these coordinate arrays are clipped so that we are in the available data
# for example, point (-1,3,7) will use the point (0,3,7) instead
m_coords = m_coords.clip(0, M.shape[0]-1)
n_coords = n_coords.clip(0, M.shape[1]-1)
o_coords = o_coords.clip(0, M.shape[2]-1)

# get the 4x4x4 cube:
A = M[m_coords, n_coords, o_coords]

# interpolate along the first axis (3D to 2D)
B = f(A[0], A[1], A[2], A[3], m_f)

# interpolate along the second axis (2D to 1D)
C = f(B[0], B[1], B[2], B[3], n_f)

# interpolate along the third axis (1D to scalar)
D = f(C[0], C[1], C[2], C[3], o_f)

return D


现在,该功能使得可以使用任何插值方法,该方法可以对四个相邻点逐个轴进行。现在显示立方,线性和“最近”。



但这有效吗?让我们通过生成一些随机数据和与数据相交的平面来对其进行测试。 (后者异常困难。)

# random data
M = np.random.random((10, 12, 16))
m0,m1 = -10.,10.
n0,n1 = -7.,5.
o0,o1 = -4.,5.

# create a grid (grid from -15..15,-15..15 in n-m plane)
gr = mgrid[-15:15.01:.1, -15:15.01:.1]

# create two perpendicular unit vectors (forming a plane)
# vn: normal vector of the plane
# vp: some vector, whose projection determines one direction
# v0: unit vector on the plane (perpendicular to vn and vp)
# v1: unit vector on the plane (perpendicular to vn and v0)
vn = np.array([-.2, .3, 1])
vp = np.array([0, -1, 0])
v1 = np.cross(vn, vp)
v2 = np.cross(vn, v1)
v1 /= numpy.linalg.norm(v1)
v2 /= numpy.linalg.norm(v2)

# the grid and the unit vectors define the 3d points on the plane
gr3d = gr[0][:,:,None] * v1 + gr[1][:,:,None] * v2

# now we can fetch the points at grid points, just must flatten and reshape back
res = [ tri_interp(M, interp_catmull_rom, p[0], p[1], p[2], m0,m1,n0,n1,o0,o1) for p in gr3d.reshape(-1,3) ]
res = np.array(res).reshape((gr3d.shape[0], gr3d.shape[1]))


因此,我们有一架穿过数据块的飞机。现在,通过选择 interp_linearinterp_nearestinterp_catmull_rom作为一维插值函数,我们实际上可以看到三种不同插值方法之间的差异。下面的所有图片在尺寸和颜色上都具有相同的比例。

1.最近:获取网格中最近的已知点(quick'n'dirty)



2.线性:取一个2x2x2立方体,并使用线性插值法计算结果;结果始终在现有数据点之间,不需要边缘欺骗。但是,第一微分不平滑,这在图像中可以看到,有时在信号处理中是个问题。



3.立方:取一个4x4x4的立方,然后使用立方插值来计算结果。结果具有连续的一阶微分,即结果在所有方向上都是“平滑的”。由于样品量大,可能会有边缘效应。 (插值算法会重复最边缘的体素,但它们可能会被镜像,取为某个常数等)。





如果我们在块中采用线段,则可以更好地看出这些插值方法的本质差异:

 # line segment from (-12,-12,-12) to (12,12,12)
v = np.array([np.linspace(-12,12,1000)]*3).T

res = [ tri_interp(M, interp_catmull_rom, p[0], p[1], p[2], m0,m1,n0,n1,o0,o1) for p in v ]


对于不同的插值函数,线段上具有以下值:



线性插值中“最近”和扭结的块状清晰可见。 “线性”插值未给出线性结果可能看起来令人惊讶。通过考虑一个带有角值为1,0,1,0的正方形的2D插值(即,对角相对的角具有相同的值)可以最好地理解。没有办法进行插值,以便所有部分都是线性的。



此处显示的算法是一种非常琐碎的算法,可以通过多种技巧使其变得更快。特别是如果在一个单元中需要多个点,则基于矩阵的方法要快得多。不幸的是,还没有非常快的方法,并且基于矩阵的方法解释起来更加复杂。

关于python - 基于网格的多元三次插值,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/24531182/

24 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com