在 x、y 和 z 中具有不同间隔的定期采样 3D 数据的快速插值
Posted
技术标签:
【中文标题】在 x、y 和 z 中具有不同间隔的定期采样 3D 数据的快速插值【英文标题】：Fast interpolation of regularly sampled 3D data with different intervals in x,y, and z 【发布时间】：2013-04-19 12:56:22 【问题描述】：我有一些体积成像数据，其中包含在 x、y、z 中的规则网格上采样的值，但具有非立方体素形状（z 中相邻点之间的空间大于 x、y 中的空间）。我最终希望能够在穿过体积的任意 2D 平面上插值，如下所示：
我知道scipy.ndimage.map_coordinates，但在我的情况下使用它并不那么简单，因为它隐含地假设输入数组中元素的间距在各个维度上是相等的。我可以首先根据最小体素维度对输入数组进行重新采样（这样我的所有体素都将成为立方体），然后使用map_coordinates 在我的平面上进行插值，但插值我的数据两次。
我还知道scipy 具有用于不规则间隔 ND 数据的各种插值器（LinearNDInterpolator、NearestNDInterpolator 等），但这些对于我的目的而言非常缓慢且占用大量内存。鉴于我知道值 在每个维度内有规律地间隔，插入数据的最佳方法是什么？
【问题讨论】：
【参考方案1】：您可以将map_coordinates 与一点代数结合使用。假设您的网格间距为dx、dy 和dz。我们需要将这些真实世界坐标映射到数组索引坐标，所以让我们定义三个新变量：
xx = x / dx
yy = y / dy
zz = z / dz
map_coordinates 的数组索引 输入是一个形状为(d, ...) 的数组，其中d 是原始数据的维数。如果你定义一个数组如：
scaling = np.array([dx, dy, dz])
您可以通过除以scaling 使用一点广播魔法将您的真实世界坐标转换为数组索引坐标：
idx = coords / scaling[(slice(None),) + (None,)*(coords.ndim-1)]
举个例子：
dx, dy, dz = 1, 1, 2
scaling = np.array([dx, dy, dz])
data = np.random.rand(10, 15, 5)
假设我们要沿平面 2*y - z = 0 插入值。我们取两个垂直于平面法线向量的向量：
u = np.array([1, 0 ,0])
v = np.array([0, 1, 2])
并得到我们想要插值的坐标：
coords = (u[:, None, None] * np.linspace(0, 9, 10)[None, :, None] +
v[:, None, None] * np.linspace(0, 2.5, 10)[None, None, :])
我们将它们转换为数组索引坐标并使用map_coordinates进行插值：
idx = coords / scaling[(slice(None),) + (None,)*(coords.ndim-1)]
new_data = ndi.map_coordinates(data, idx)
最后一个数组的形状为(10, 10)，并且在[u_idx, v_idx] 位置具有对应于坐标coords[:, u_idx, v_idx] 的值。
您可以基于这个想法来处理坐标不从零开始的插值，方法是在缩放之前添加一个偏移量。
【讨论】：
这正是我所需要的。干杯，詹姆！【参考方案2】：这是一个简单的类Intergrid
将非均匀映射/缩放到均匀网格，
然后map_coordinates.
在4d test case 上，它以每个查询点大约 1 微秒的速度运行。
pip install [--user] intergrid 应该可以在 python2 或 python3 中工作（2020 年 2 月）；见intergrid on PyPi。
""" interpolate data given on an Nd rectangular grid, uniform or non-uniform.
Purpose: extend the fast N-dimensional interpolator
`scipy.ndimage.map_coordinates` to non-uniform grids, using `np.interp`.
Background: please look at
http://en.wikipedia.org/wiki/Bilinear_interpolation
https://***.com/questions/6238250/multivariate-spline-interpolation-in-python-scipy
http://docs.scipy.org/doc/scipy-dev/reference/generated/scipy.ndimage.interpolation.map_coordinates.html
Example
-------
Say we have rainfall on a 4 x 5 grid of rectangles, lat 52 .. 55 x lon -10 .. -6,
and want to interpolate (estimate) rainfall at 1000 query points
in between the grid points.
# define the grid --
griddata = np.loadtxt(...) # griddata.shape == (4, 5)
lo = np.array([ 52, -10 ]) # lowest lat, lowest lon
hi = np.array([ 55, -6 ]) # highest lat, highest lon
# set up an interpolator function "interfunc()" with class Intergrid --
interfunc = Intergrid( griddata, lo=lo, hi=hi )
# generate 1000 random query points, lo <= [lat, lon] <= hi --
query_points = lo + np.random.uniform( size=(1000, 2) ) * (hi - lo)
# get rainfall at the 1000 query points --
query_values = interfunc( query_points ) # -> 1000 values
What this does:
for each [lat, lon] in query_points:
1) find the square of griddata it's in,
e.g. [52.5, -8.1] -> [0, 3] [0, 4] [1, 4] [1, 3]
2) do bilinear (multilinear) interpolation in that square,
using `scipy.ndimage.map_coordinates` .
Check:
interfunc( lo ) -> griddata[0, 0],
interfunc( hi ) -> griddata[-1, -1] i.e. griddata[3, 4]
Parameters
----------
griddata: numpy array_like, 2d 3d 4d ...
lo, hi: user coordinates of the corners of griddata, 1d array-like, lo < hi
maps: a list of `dim` descriptors of piecewise-linear or nonlinear maps,
e.g. [[50, 52, 62, 63], None] # uniformize lat, linear lon
copy: make a copy of query_points, default True;
copy=False overwrites query_points, runs in less memory
verbose: default 1: print a 1-line summary for each call, with run time
order=1: see `map_coordinates`
prefilter: 0 or False, the default: smoothing B-spline
1 or True: exact-fit interpolating spline (IIR, not C-R)
1/3: Mitchell-Netravali spline, 1/3 B + 2/3 fit
(prefilter is only for order > 1, since order = 1 interpolates)
Non-uniform rectangular grids
-----------------------------
What if our griddata above is at non-uniformly-spaced latitudes,
say [50, 52, 62, 63] ? `Intergrid` can "uniformize" these
before interpolation, like this:
lo = np.array([ 50, -10 ])
hi = np.array([ 63, -6 ])
maps = [[50, 52, 62, 63], None] # uniformize lat, linear lon
interfunc = Intergrid( griddata, lo=lo, hi=hi, maps=maps )
This will map (transform, stretch, warp) the lats in query_points column 0
to array coordinates in the range 0 .. 3, using `np.interp` to do
piecewise-linear (PWL) mapping:
50 51 52 53 54 55 56 57 58 59 60 61 62 63 # lo[0] .. hi[0]
0 .5 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 3
`maps[1] None` says to map the lons in query_points column 1 linearly:
-10 -9 -8 -7 -6 # lo[1] .. hi[1]
0 1 2 3 4
More doc: https://denis-bz.github.com/docs/intergrid.html
"""
# split class Gridmap ?
from __future__ import division
from time import time
# warnings
import numpy as np
from scipy.ndimage import map_coordinates, spline_filter
__version__ = "2014-01-15 jan denis" # 15jan: fix bug in linear scaling
__author_email__ = "denis-bz-py@t-online.de" # comments welcome, testcases most welcome
#...............................................................................
class Intergrid:
__doc__ = globals()["__doc__"]
def __init__( self, griddata, lo, hi, maps=[], copy=True, verbose=1,
order=1, prefilter=False ):
griddata = np.asanyarray( griddata )
dim = griddata.ndim # - (griddata.shape[-1] == 1) # ??
assert dim >= 2, griddata.shape
self.dim = dim
if np.isscalar(lo):
lo *= np.ones(dim)
if np.isscalar(hi):
hi *= np.ones(dim)
self.loclip = lo = np.asarray_chkfinite( lo ).copy()
self.hiclip = hi = np.asarray_chkfinite( hi ).copy()
assert lo.shape == (dim,), lo.shape
assert hi.shape == (dim,), hi.shape
self.copy = copy
self.verbose = verbose
self.order = order
if order > 1 and 0 < prefilter < 1: # 1/3: Mitchell-Netravali = 1/3 B + 2/3 fit
exactfit = spline_filter( griddata ) # see Unser
griddata += prefilter * (exactfit - griddata)
prefilter = False
self.griddata = griddata
self.prefilter = (prefilter == True)
self.maps = maps
self.nmap = 0
if len(maps) > 0:
assert len(maps) == dim, "maps must have len %d, not %d" % (
dim, len(maps))
# linear maps (map None): Xcol -= lo *= scale -> [0, n-1]
# nonlinear: np.interp e.g. [50 52 62 63] -> [0 1 2 3]
self._lo = np.zeros(dim)
self._scale = np.ones(dim)
for j, (map, n, l, h) in enumerate( zip( maps, griddata.shape, lo, hi )):
## print "test: j map n l h:", j, map, n, l, h
if map is None or callable(map):
self._lo[j] = l
if h > l:
self._scale[j] = (n - 1) / (h - l) # _map lo -> 0, hi -> n - 1
else:
self._scale[j] = 0 # h <= l: X[:,j] -> 0
continue
self.maps[j] = map = np.asanyarray(map)
self.nmap += 1
assert len(map) == n, "maps[%d] must have len %d, not %d" % (
j, n, len(map) )
mlo, mhi = map.min(), map.max()
if not (l <= mlo <= mhi <= h):
print "Warning: Intergrid maps[%d] min %.3g max %.3g " \
"are outside lo %.3g hi %.3g" % (
j, mlo, mhi, l, h )
#...............................................................................
def _map_to_uniform_grid( self, X ):
""" clip, map X linear / nonlinear inplace """
np.clip( X, self.loclip, self.hiclip, out=X )
# X nonlinear maps inplace --
for j, map in enumerate(self.maps):
if map is None:
continue
if callable(map):
X[:,j] = map( X[:,j] ) # clip again ?
else:
# PWL e.g. [50 52 62 63] -> [0 1 2 3] --
X[:,j] = np.interp( X[:,j], map, np.arange(len(map)) )
# linear map the rest, inplace (nonlinear _lo 0, _scale 1: noop)
if self.nmap < self.dim:
X -= self._lo
X *= self._scale # (griddata.shape - 1) / (hi - lo)
## print "test: _map_to_uniform_grid", X.T
#...............................................................................
def __call__( self, X, out=None ):
""" query_values = Intergrid(...) ( query_points npt x dim )
"""
X = np.asanyarray(X)
assert X.shape[-1] == self.dim, ("the query array must have %d columns, "
"but its shape is %s" % (self.dim, X.shape) )
Xdim = X.ndim
if Xdim == 1:
X = np.asarray([X]) # in a single point -> out scalar
if self.copy:
X = X.copy()
assert X.ndim == 2, X.shape
npt = X.shape[0]
if out is None:
out = np.empty( npt, dtype=self.griddata.dtype )
t0 = time()
self._map_to_uniform_grid( X ) # X inplace
#...............................................................................
map_coordinates( self.griddata, X.T,
order=self.order, prefilter=self.prefilter,
mode="nearest", # outside -> edge
# test: mode="constant", cval=np.NaN,
output=out )
if self.verbose:
print "Intergrid: %.3g msec %d points in a %s grid %d maps order %d" % (
(time() - t0) * 1000, npt, self.griddata.shape, self.nmap, self.order )
return out if Xdim == 2 else out[0]
at = __call__
# end intergrid.py
【讨论】：
【参考方案3】：我创建了 regulargrid 包（https://pypi.python.org/pypi/regulargrid/，来源https://github.com/JohannesBuchner/regulargrid）
它通过非常快速的 scipy.ndimage.map_coordinates 为任意坐标比例提供对 n 维笛卡尔网格（根据需要）的支持。
另请参阅此答案：Fast interpolation of grid data
【讨论】：
以上是关于在 x、y 和 z 中具有不同间隔的定期采样 3D 数据的快速插值的主要内容，如果未能解决你的问题，请参考以下文章