数据集的高效拼接

Question

我有多个测量数据集，我想将它们合并为一个数据集。虽然我有一个可行的解决方案，但效率非常低，我很乐意提供一些有关如何改进它的提示。

将测量视为一个对象的多个高度图，我想将其合并为一个高度图。我的测量结果并不完美，可能会有一些倾斜和高度偏移。让我们假设（现在）我们知道 x-y 位置非常准确。这是一个例子：

import numpy as np
import matplotlib.pyplot as plt

def height_profile(x, y):
    radius = 100
    return np.sqrt(radius**2-x**2-y**2)-radius

np.random.seed(123)

datasets = {}

# DATASET 1
x = np.arange(-8, 2.01, 0.1)
y = np.arange(-3, 7.01, 0.1)

xx, yy = np.meshgrid(x, y)
# height is the actual profile + noise
zz = height_profile(xx, yy) + np.random.randn(*xx.shape)*0.001

datasets[1] = [xx, yy, zz]

plt.figure()
plt.pcolormesh(*datasets[1])
plt.colorbar()

# DATASET 2
x = np.arange(-2, 8.01, 0.1)
y = np.arange(-3, 7.01, 0.1)

xx, yy = np.meshgrid(x, y)
# height is the actual profile + noise + random offset + random tilt
zz = height_profile(xx, yy) + np.random.randn(*xx.shape)*0.001 + np.random.rand() + np.random.rand()*xx*0.1 + np.random.rand()*yy*0.1

datasets[2] = [xx, yy, zz]

plt.figure()
plt.pcolormesh(*datasets[2])
plt.colorbar()

# DATASET 3
x = np.arange(-5, 5.01, 0.1)
y = np.arange(-7, 3.01, 0.1)

xx, yy = np.meshgrid(x, y)
# height is the actual profile + noise + random offset + random tilt
zz = height_profile(xx, yy) + np.random.randn(*xx.shape)*0.001 + np.random.rand() + np.random.rand()*xx*0.1 + np.random.rand()*yy*0.1

datasets[3] = [xx, yy, zz]

plt.figure()
plt.pcolormesh(*datasets[3])
plt.colorbar()

为了组合三个（或更多）数据集，我有以下策略：找到数据集之间的重叠，计算重叠区域（residual_overlap）中数据集之间的总高度差并尝试使用 lmfit 最小化高度差（剩余）。要在数据集上应用 t运行sformations（倾斜、偏移等），我有一个专门的功能。

from lmfit import minimize, Parameters
from copy import deepcopy
from itertools import combinations
from scipy.interpolate import griddata

def data_transformation(dataset, idx, params):
    dataset = deepcopy(dataset)
    
    if 'x_offset_{}'.format(idx) in params:
        x_offset = params['x_offset_{}'.format(idx)].value
    else:
        x_offset = 0

    if 'y_offset_{}'.format(idx) in params:
        y_offset = params['y_offset_{}'.format(idx)].value
    else:
        y_offset = 0
    
    if 'tilt_x_{}'.format(idx) in params:
        x_tilt = params['tilt_x_{}'.format(idx)].value
    else:
        x_tilt = 0

    if 'tilt_y_{}'.format(idx) in params:
        y_tilt = params['tilt_y_{}'.format(idx)].value
    else:
        y_tilt = 0

    if 'piston_{}'.format(idx) in params:
        piston = params['piston_{}'.format(idx)].value
    else:
        piston = 0

    _x = dataset[0] - np.mean(dataset[0])
    _y = dataset[1] - np.mean(dataset[1])

    dataset[0] = dataset[0] + x_offset
    dataset[1] = dataset[1] + y_offset
    dataset[2] = dataset[2] + 2 * (x_tilt * _x + y_tilt * _y) + piston

    return dataset

def residual_overlap(dataset_0, dataset_1):
    xy_0 = np.stack((dataset_0[0].flatten(), dataset_0[1].flatten()), axis=1)
    xy_1 = np.stack((dataset_1[0].flatten(), dataset_1[1].flatten()), axis=1)
    difference = griddata(xy_0, dataset_0[2].flatten(), xy_1) - \
                 dataset_1[2].flatten()

    return difference

def residual(params, datasets):
    datasets = deepcopy(datasets)

    for idx in datasets:
        datasets[idx] = data_transformation(
            datasets[idx], idx, params)

    residuals = []

    for combination in combinations(list(datasets), 2):
        residuals.append(residual_overlap(
            datasets[combination[0]], datasets[combination[1]]))

    residuals = np.concatenate(residuals)
    residuals[np.isnan(residuals)] = 0

    return residuals

def minimize_datasets(params, datasets, **minimizer_kw):
    minimize_fnc = lambda *args, **kwargs: residual(*args, **kwargs)

    datasets = deepcopy(datasets)

    min_result = minimize(minimize_fnc, params,
                          args=(datasets, ), **minimizer_kw)

    return min_result

我运行这样的“拼接”：

params = Parameters()
params.add('tilt_x_2', 0)
params.add('tilt_y_2', 0)
params.add('piston_2', 0)
params.add('tilt_x_3', 0)
params.add('tilt_y_3', 0)
params.add('piston_3', 0)

fit_result = minimize_datasets(params, datasets)

plt.figure()
plt.pcolormesh(*data_transformation(datasets[1], 1, fit_result.params), alpha=0.3, vmin=-0.5, vmax=0)
plt.pcolormesh(*data_transformation(datasets[2], 2, fit_result.params), alpha=0.3, vmin=-0.5, vmax=0)
plt.pcolormesh(*data_transformation(datasets[3], 3, fit_result.params), alpha=0.3, vmin=-0.5, vmax=0)
plt.colorbar()

如您所见，它确实有效，但在我的计算机上拼接这些小数据集大约需要一分钟。实际上我有更多更大的数据集。

您是否找到提高拼接性能的方法？

编辑： 按照建议，我运行一个分析器，它显示 99.5% 的时间花在了 griddata 函数上。那个用于将数据点从 dataset_0 插入到 dataset_1 的位置。如果我将方法切换到“最近”，执行时间会下降到大约一秒，但不会发生插值。有机会提高插值速度吗？

Answer 1

浏览代码，除了你一遍又一遍地运行 deepcopy() 之外，我真的看不到任何可以改进的地方。

但是，我建议您这样做 profiling。如果您使用 pycharm，则可以使用 clock/run 符号进行分析。

我相信其他 IDE 也有这样的能力。这样你就可以找出哪个函数花费的时间最多。

整个图表：

当我放大几个功能时（我正在显示 google 云功能）：

你可以看到他们被调用了多少次，他们花了多长时间等

长话短说，你需要一个分析器！

数据集的高效拼接

Efficient stitching of datasets

python

optimization

lmfit