python - 构建一个 2D 数组，表示由 3 个点和所需的数组大小定义的 3D 平面（存储其 Z 值）

Question

短的

我需要找到一种优化的方法来构建一个表示 3D 平面（即存储其 Z 值）的 2D numpy 数组，给定三个点（以数学方式定义平面）和所需的数组大小。

细节

给定三个点p0、p1、p2和一个size变量：

import numpy as np
import pylab as plt

p0 = [0.2, -0.4, -0.2]
p1 = [-0.6, 0.1, -0.8]
p2 = [-0.1, 0.3, -0.6]
size = (12000, 17000)

以下函数显示了如何从这 3 个点构建平面的想法：

def plane(p0, p1, p2, size):
    p0 = np.array(p0)
    p1 = np.array(p1)
    p2 = np.array(p2)
    ux, uy, uz = [p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]]
    vx, vy, vz = [p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]]
    uv = [uy * vz - uz * vy,
          uz * vx - ux * vz,
          ux * vy - uy * vx]
    norm = np.array(uv)    
    d = -np.array(p0).dot(norm)    
    xx, yy = np.meshgrid(range(size[1]), range(size[0]))    
    z = np.array((-norm[0] * xx - norm[1] * yy - d) * 1. / norm[2])
    return(z)

结果z如下所示：

这个实现会在几秒钟内吃掉我的记忆：

问题

您将如何优化它以尽可能低的内存使用达到完全相同的结果（快速也不是可选的......）？

注意：size可能在两个维度上都大 10 倍 -> 这显然会使我的 Python 在 2016 年主流笔记本电脑上崩溃。

（使用最常见的科学图书馆完全没问题）

我在 Ubuntu linux 18.04.5 上使用 Python 3.6.9（默认，2021 年 1 月 26 日，15:33:00）[GCC 8.4.0]。

标签： pythonnumpyoptimizationmemoryplane

Answer 1

在 的计算之前d，您的数组很小。之后的一切都是为操作的每个部分创建一个单独的 1.5GiB (12k * 17k * 8bytes/float) 临时数组。您可以完全摆脱这些阵列中的大部分，并在原地完成其余的工作。例如：

xx, yy = np.meshgrid(...)

您真的不需要 meshgrid: 广播np.arange(size[0])[:, None]并将np.arange(size[1])制作您需要的所有阵列。对于这部分数据，这会将您的大小从 ~3GiB 减少到 ~230KiB：

xx = np.arange(size[1], dtype=float)
yy = np.arange(size[0], dtype=float)[:, None]

该表达式z = ...中大约有七个临时数组，每个数组大小为 1.5GiB：

1. t1 = -norm[0] * xx
2. t2 = norm[1] * yy
3. t3 = t1 - t2
4. t4 = t3 - d
5. t5 = t4 * 1
6. t6 = t5 / norm[2]
7. t7 = np.array(t6)

其中一些数组可能在计算过程中被垃圾收集，但它们都需要在某个时候被分配。您可以使用广播和就地运算符或out适当ufuncs的参数来摆脱所有这些：

所以：

d /= norm[2]
norm /= norm[2]
xx *= -norm[0]
yy *= -norm[1]
yy -= d
z = xx + yy

即使在代码的紧凑部分，你也在做一堆不必要的操作。例如，你打np.array(p0)了两次电话。这是您功能的更好版本：

def plane2(p0, p1, p2, size):
    u = np.subtract(p1, p0)
    v = np.subtract(p2, p0)
    uxv = np.cross(u, v)
    d = -np.dot(p0, uxv)
    d /= uxv[2]
    uxv /= uxv[2]
    xx = np.arange(size[1], dtype=float)
    xx *= -uxv[0]
    yy = np.arange(size[0], dtype=float)[:, None]
    yy *= -uxv[1]
    yy -= d
    return xx + yy

对于比您建议的尺寸小 10 倍的尺寸：

%timeit plane(p0, p1, p2, size)
70.1 ms ± 329 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit plane2(p0, p1, p2, size)
8.29 ms ± 302 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

当我将大小增加到您问题中的大小时，您的实现会崩溃，而我的实现会在合理的时间内完成：

%timeit plane(p0, p1, p2, size)
...
MemoryError: Unable to allocate 1.52 GiB for an array with shape (12000, 17000) and data type float64

%timeit plane2(p0, p1, p2, size)
819 ms ± 37.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

根据您拥有多少 RAM 以及 python 进程中的可用内存设置，您可能无法将数组大小增加 100 倍至 152GiB，无论您执行操作的效率如何。