用于在 3D 点阵中查找 2D 点阵的拓扑匹配算法

Question

我有一个 3D 晶格，其单位单元（即最小重复单位）为 16 个点。由于它是 3D 晶格，因此它在所有 3 个维度（x、y、z）上都是周期性的。详细来说，单位单元如下所示：\n

\n

3 个晶格向量 a ₁、a ₂、a ₃为

\n

bulk_vecs = [[14.56578026795, 0.0,            0.0          ], # a1\n             [0.0,            8.919682340494, 0.0          ], # a2\n             [7.282890133975, 2.973227446831, 4.20477857933]] # a3\n

\n

16个点的笛卡尔坐标为

\n

bulk_coords = [[ 0.00000000,  0.00000000,  0.00000000],\n               [ 10.9243352,  6.34420137,  2.66488775],\n               [ 18.2072253,  6.34420137,  2.66488775],\n               [ 12.7450577,  3.17210068,  1.33244387],\n               [ 10.6807662,  8.91968234,  0.42187384],\n               [ 16.1429338,  3.37097392,  1.19181926],\n               [ 19.7843789,  9.31742882,  3.29420855],\n               [ 9.34718165,  2.97322745,  1.47306849],\n               [ 14.8093492,  6.34420137,  2.24301390],\n               [ 9.10361267,  9.11855558,  3.43483316],\n               [ 3.88501404,  0.39774648,  0.14062462],\n               [ 7.03932116,  5.94645489,  2.52426313],\n               [ 12.9886267,  8.91968234,  3.57545778],\n               [ 7.28289013,  0.00000000,  0.00000000],\n               [ 5.46216760,  3.17210068,  1.33244387],\n               [ 16.3865028,  9.11855558,  3.43483316]]\n

\n

现在我有一个参考 2D 晶格，其中包含单位单元中的 32 个点。请注意，这里的“2D”意味着它仅在 x 和 y 维度上是周期性的，但在 z 方向上仍然可以具有厚度。它的晶胞看起来像这样（这只是一个简单的例子，晶胞可以是任何形状，不一定是立方体）：

\n

\n

3 个晶格向量由下式给出

\n

ref_vecs = [[8.1968234, 0.,          0.          ], # x-direction (periodic)\n            [0.,        13.38535656, 0.          ], # y-direction (periodic)\n            [0.,        0.,          7.7280392142]] # z-thickness (non-periodic)\n

\n

32个点的笛卡尔坐标为

\n

ref_coords = [[ 5.65852755,  9.84826406,  0.10024035],\n              [ 5.66769587,  3.14583318,  0.03049278],\n              [ 5.84383908,  6.16622816,  0.33687635],\n              [ 3.06154746,  11.5036515,  0.89497370],\n              [ 3.04533245,  4.74684988,  0.80482405],\n              [ 2.81714593,  7.85388222,  0.95654678],\n              [ 3.09409158,  1.66514480,  1.21902173],\n              [-0.03484920,  9.98872688,  1.91238251],\n              [ 0.07586672,  6.55689489,  2.02252838],\n              [ 0.15624468,  13.1218508,  2.07583316],\n              [ 0.30441248,  2.82022739,  2.21046382],\n              [ 5.46226670,  11.2187312,  2.95396399],\n              [ 5.83993074,  5.03326172,  3.11251461],\n              [ 5.40963696,  8.13117939,  3.23409626],\n              [ 5.45612978,  1.47718658,  3.32634157],\n              [ 2.72591886,  6.68506441,  3.86756751],\n              [ 2.91602855,  3.06649570,  3.90688764],\n              [ 2.95946887,  9.79329729,  3.96116732],\n              [ 3.11916042,  12.9792231,  4.17677610],\n              [ 0.35863087,  4.72533071,  4.75876360],\n              [ 0.31594175,  11.4937246,  4.67355664],\n              [ 0.03089551,  1.37091216,  4.91348486],\n              [ 0.39020723,  8.35223658,  5.11836201],\n              [ 5.45925351,  3.34155067,  5.87036690],\n              [ 5.62981527,  13.2212649,  5.93479016],\n              [ 5.64259931,  6.46196620,  5.93435713],\n              [ 5.83825398,  9.54653517,  6.08472919],\n              [ 2.80592452,  4.61964086,  6.76956275],\n              [ 3.15087591,  11.6754589,  6.98585963],\n              [ 2.68542186,  1.40293444,  7.10791606],\n              [ 2.73378423,  8.13529985,  7.12898617],\n              [ 0.03114729, -0.01057378,  7.78083813]]\n

\n

该参考 2D 晶格表示类似于从 3D 晶格中提取的“厚层”的东西，并且其 x 和 y 方向上的 2 个晶格向量平行于a 1_和a 2 _（即平行于 xy 平面）。现在我想在当前 3D 晶格内找到与该参考 2D 晶格最匹配的并行 2D 晶格。相似性应取决于两个二维晶格晶胞​​的拓扑结构。例如，参考晶胞可以被认为是对最佳匹配的二维晶格晶胞​​的分数坐标和晶格参数应用小扰动。下面是从 3D 晶格提取 2D 晶格的示意图：\n

\n

在 3D 晶格空间中搜索所有可能的并行 2D 晶格单元的空间是很困难的，特别是考虑到我可能需要扩展 3D 晶格单元以在整个空间中搜索（但我应该扩展多少？）。我想知道是否有任何算法可以在整个3D晶格空间中相对较快地找到最佳拓扑匹配的2D晶格单元\xef\xbc\x9f具体来说，我想获得晶格向量和32最佳匹配的二维晶格晶胞​​的笛卡尔坐标。欢迎任何建议甚至想法。

\n

编辑：

\n

好的，我会尝试以最简单的方式解释这个问题。正如您可能已经猜到的那样，我正在研究晶体学。3D 晶格表示体晶体结构（即所有 3 个维度上的周期性原子排列），2D 晶格表示表面结构（即 x 和 y 方向上的周期性原子排列，但在 z 方向上具有非周期性厚原子层） ）。显然，可以通过切割块体结构来获得表面结构。

\n

让我们使用一个简单的 Si 体结构（从MaterialsProject下载）作为示例，它有一个立方晶胞，有 3 个晶格向量

\n

bulk_vecs = [[5.44370237, 0.0,        0.0       ],\n             [0.0,        5.44370237, 0.0       ],\n             [0.0,        0.0,        5.44370237]]\n

\n

和原子坐标

\n

bulk_coords = [[ 1.27013575,  0.03307657, -0.02898786],\n               [ 0.02039057,  0.02332397,  0.79145084],\n               [ 0.0137774,   1.18012325,  1.85879242],\n               [ 1.22267376,  1.27420799,  2.72875596],\n               [ 1.20645874, -0.0488289,   3.52879259],\n               [ 0.08487146, -0.06430577,  4.44147333]]\n

\n

金刚石具有与硅相同类型的晶体结构，但晶格参数和原子坐标略有不同。现在假设我有一个参考钻石 (001) 表面结构（米勒指数(001) 基本上意味着表面平面平行于 xy 平面）。它有一个 6 原子晶胞，3 个晶格向量为

\n

ref_vecs = [[2.5178270133578393, 0.0,                0.0        ],\n            [0.0,                2.5178270133578393, 0.0        ],\n            [0.0,                0.0,                5.341117665]]\n

\n

和原子坐标

\n

ref_coords = [[ 1.27013575,  0.03307657, -0.02898786],\n              [ 0.02039057,  0.02332397,  0.79145084],\n              [ 0.0137774,   1.18012325,  1.85879242],\n              [ 1.22267376,  1.27420799,  2.72875596],\n              [ 1.20645874, -0.0488289,   3.52879259],\n              [ 0.08487146, -0.06430577,  4.44147333]]\n

\n

由于金刚石具有与 Si 相同类型的晶体结构，因此应该能够从参考金刚石 (001) 表面晶胞中找到 Si(001) 表面晶胞。据我所知，一个简单的拓扑匹配算法可能是在大块 Si 空间中找到一个 Si(001) 表面晶胞，以便与参考菱形 (001) 表面晶胞的分数坐标差之和被最小化。然而，我不知道如何在这么大的空间中进行搜索，特别是当参考表面晶胞很大时，可能需要将体晶胞扩大很多，这意味着搜索空间变得更大。

\n

出于基准测试的目的，我将为最匹配的 Si(001) 表面晶胞提供预期的解决方案。它的晶格向量应该接近

\n

sol_vecs = [[3.8492788605882797, 0.0,                0.0        ],\n            [0.0,                3.8492788605882797, 0.0        ],\n            [0.0,                0.0,                8.165553555]]\n

\n

和原子坐标接近

\n

sol_coords = [[ 1.92463943, 0.0,        0.0       ],\n              [ 0.0,        0.0,        1.36092562],\n              [ 0.0,        1.92463943, 2.72185121],\n              [ 1.92463943, 1.92463943, 4.08277680],\n              [ 1.92463943, 0.0,        5.44370240],\n              [ 0.0,        0.0,        6.80462799]]\n

\n

这是晶体学方面这个问题的说明（由于某种原因无法上传到 SO）。如图所示，我提供的 Si(001) 表面晶胞溶液与参考钻石(001) 表面晶胞非常匹配。

\n

提示：在倒易空间上工作可能比在真实空间上工作容易得多。

\n

最终编辑：

\n

好的，让我以最清晰和简单的方式重新表述这个问题：如何在整个 3D 晶格空间中找到一个单元，其中它包含与参考 2D 晶格单位单元相同数量的点，并且在两个细胞的分数坐标（即按细胞大小归一化的笛卡尔坐标）？我们拥有的信息是 (1) 3D 晶格晶胞 (2) 参考 2D 晶格晶胞 (3) 2D 晶格平行于 3D 晶格的 xy 平面

\n

Answer 1

介绍

\n

让我稍微重新表述一下这个问题：我们有一个用点填充空间的 3D 图案。我们有一个参考平面，基本上可以在任何地方切割这个空间，但是我们只知道空间中的点在平面上创建的图案，没有其他信息。\n找到与参考平面最匹配的平面。\n更糟糕的是，图案不是完全二维的，它更类似于薄片。

\n

对我们来说幸运的是，点数相当有限，因此我们可以采用一些暴力方法。然而，我们不能仅仅计算所有可能的二维晶格并找到最匹配的一个，组合数学会杀了我们。但我们可以做类似的事情，只是更复杂一些。

\n

设置游乐场：

\n

首先，让我们生成我们将要工作的空间。为了简单起见，让我们沿着所有轴延伸 3D 晶格，以生成一个足够大的空间，以适应 2D 晶格的任何给定方向和位移。\n在我们的玩具示例中， 3D 晶格向量为我们提供了一个立方体，而 2D 晶格向量适合该立方体。因此，如果我们在每一侧添加一些填充，这将覆盖所有可能的旋转。现在，如果参考 2D 晶格很大或者 3D 晶格倾斜，这可能还不够，因此请记住\xe2\x80\x94，它可能需要一些数学运算并以创建足够大空间的方式堆叠晶格适合参考二维晶格的所有可能的旋转和平移。
\n如果我们错过了中心 3D 格子，我们总是可以通过声明与我们有最大重叠的格子来快速回到它，这是新的中心部分。

\n

希望不会太多而使我们的后续工作无效。

\n

笔记：

\n

如果一个格子有 32 个点，那么我们最多可以使用 300 个格子来创建我们的空间。到目前为止，匹配算法最糟糕的部分是二次的，因此 10K * 10K * 2d_lattice_points 计算应该是可行的。300 个格子给我们提供了类似 6x6x6 格子立方体的东西。另外，由于我们担心“扁平且尖”的情况，因此我们正在寻找一些矩形（例如，如果晶格是 3x3x10，那么我们可以进行 10x10x3 堆叠来补偿）。

\n

笔记：

\n

如果 3D 晶格倾斜，我建议剪掉倾斜的部分并在其前面创建一个矩形。由于周期性，这应该没问题。

\n

import numpy as np\n\nbulk_vecs = np.array([[5.44370237, 0.0,        0.0       ],\n                      [0.0,        5.44370237, 0.0       ],\n                      [0.0,        0.0,        5.44370237]])\n\nbulk_coords = np.array([[ 1.27013575,  0.03307657, -0.02898786],\n                       [ 0.02039057,  0.02332397,  0.79145084],\n                       [ 0.0137774,   1.18012325,  1.85879242],\n                       [ 1.22267376,  1.27420799,  2.72875596],\n                       [ 1.20645874, -0.0488289,   3.52879259],\n                       [ 0.08487146, -0.06430577,  4.44147333]])\n\nref_vecs = np.array([[2.5178270133578393, 0.0,                0.0        ],\n                    [0.0,                2.5178270133578393, 0.0        ],\n                    [0.0,                0.0,                5.341117665]])\n\nbulk_space = []\nfor point in [x for x in bulk_coords]:\n    for i in range(3):\n        for j in range(3):\n            for k in range(3):\n                shift = i * bulk_vecs[0] + j * bulk_vecs[1] + k * bulk_vecs[2]\n                bulk_space.append(point + shift)\n

\n

这在玩具示例中为我们提供了大约 162 个点，这很好。如果这个级别的堆叠对于主要情况来说足够了，那么这会给我们一些 432 点，这仍然可以。我的希望是进行足够的优化，以便能够在大空间中达到数千个点。（例如，3D 晶格生成的游乐场。）
\n超出这个范围可能是可能的，但很快就会变得烦人\xe2\x80\x94 正如我提到的，我将进行一些部分暴力破解。

\n

让巫毒开始吧

\n

首先，让我们去掉所有那些讨厌的小数。计算机不喜欢它们。\n为了做到这一点，让我们同意，三位小数足以满足精度（目前）

\n

bulk_space = np.array(np.array(bulk_space) * 1000, dtype='longlong')\n

\n

请注意，这可能太少或太多\xe2\x80\x94我不确定 2D 平面内发生的变形/不精确的实际水平。
\n我会尝试做一些匹配\xe2\x80\x94如果匹配太多我会提高精度。如果太少（或没有）我会减少它。
\n这是最好的离散化\xe2\x80\x94我们将连续空间（或更精细，在这种情况下我们已经没有完整的连续体）分割成固定大小的框，然后说每个点都与一个框匹配它包含在其中。
\n然后我们忘记了所有原始点并愉快地用盒子进行计算。
\n最后，如果我们的努力成功了，我们就会回到要点并尝试提高精度。

\n

笔记：

\n

从我得到的情况来看，所选的精度似乎不会导致点（立方体内的两个点）\xe2\x80\x94的任何重叠，我们正在工作的空间似乎人口密度不太高。如果两个（或更多）点重叠，我们不会太介意。然而，这意味着在最终检查中（同时提高精度）我们需要用我们的解决方案检查每一个。

\n

平面近似：

\n

由于我们不了解 2D 平面方向，因此我们将首先使用距离匹配。我们将计算 3D 空间的距离矩阵并将其存储起来。

\n

然后我们将计算 2D 晶格的距离矩阵，并取最长距离，因为由于 3D 晶格盒中的空间有限，它的重复次数最少。

\n

然后我们将取它的一个端点，并取到其他点的所有距离并对它们进行排序。

\n

笔记：

\n

我们将比较距离，我们实际上并不关心它们的值，只关心它们的比较，我们可以省略平方根并使用距离的平方。它的计算成本更便宜。

\n

现在我们将搜索 2D 和 3D 中最大距离的匹配。对于每场比赛，我们将尝试找到匹配的距离序列。

\n

笔记2：

\n

为了检查匹配是否存在，我们可以对与该点对应的距离数组进行排序，找到匹配的最大值（我们已经知道它从前一个点就存在），然后按顺序找到其余的\xe2\x80\x94如果我们找不到任何给定的匹配项，我们可以丢弃整个序列，因为它已排序。
\n请注意，我们要么需要在排序时保留索引（这有点烦人），要么在确认匹配项存在后对匹配项进行暴力搜索，以从矩阵推断索引（这需要一点时间\xe2\x80\x94，但不需要太多了，只是基于 2D 的点数乘以 3D 点数）。
\n所以我们会这么做。

\n

注3：

\n

如果由于二维晶格切割的不精确而没有找到任何匹配，我们可以以较低的精度再次尝试，直到找到它。

\n

注4：

\n

此时我们不关心找到多个匹配项。我们只是过滤/验证这些点。我们将在下一步中处理这些问题。

\n

for i in range(len(ref_space)):\n    for j in range(i, len(ref_space)):\n        diff = ref_space[i] - ref_space[j]\n        dist = np.sum(diff * diff)\n        distance_matrix_2d[i, j] = dist\n        distance_matrix_2d[j, i] = dist\n        if dist > distance_max_2d:\n            dmax_i = i\n            dmax_j = j\n            distance_max_2d = dist\n\ndistance_matrix_3d = np.zeros(shape=[len(bulk_space), len(bulk_space)], dtype='longlong')\nmatch_list_i = []\nmatch_list_j = []\npoint_match = []\ndist_min = 1000000000\ndist_max = 0\nfor i in range(len(bulk_space)):\n    for j in range(i, len(bulk_space)):\n        diff = bulk_space[i] - bulk_space[j]\n        dist = np.sum(diff * diff)\n        distance_matrix_3d[i, j] = dist\n        distance_matrix_3d[j, i] = dist\n        if distance_max_2d - dist == 0:\n            match_list_i.append(i)\n            match_list_j.append(j)\n            point_match.append(i)\n            point_match.append(j)\n\npoint_matches = []\nfor i in point_match:\n    x = []\n    for j in range(len(distance_matrix_3d)):\n        x.append(distance_matrix_3d[i, j])\n    x = np.sort(x)\n    point_matches.append((i, x))\n# Note that we're currently not using the sort\n# and I'm simply brute forcing the next step.\n# If that is an issue, this sort can be used to\n# make the next step log(n) * 2d_pts_count\n# by checking first the existence of the\n# distance via binary search in the 3D\n# sorted matrix.\n# However since we only have few 2D points, this\n# is likely unnecessary unless we're solving\n# some big case, where there might be other\n# performance bottlenecks.\n\npoint_matches_bruteforce = []\nfor i, matches in point_matches:\n    # Check all the rows defined by points from the previous\n    # step (selected by matching the longest distance)\n    row_failed = False\n    for dist2d in distance_matrix_2d[dmax_i]:\n        # For every distance in 2D matrix,\n        # try finding a matching one\n        found = False\n        for dist3d in matches:\n            # If we find it, we qualify the distance\n            if dist3d == dist2d:\n                found = True\n        #print("dist2d checked")\n        if not found:\n            # If we didn't find any match, we can discard\n            # the point. There is at least one point from\n            # the set without matching distance, e.g.,\n            # we can't match the structure.\n            row_failed = True\n            break\n    if row_failed:\n        # If any row failed (no matching 3D distance\n        # for a 2D one, we discard the point)\n        print("point ", i, "has failed to find matching set")\n    else:\n        # Otherwise, we keep it for the next, more bruteforcish round\n        print("point ", i, "is a candidate")\n        point_matches_bruteforce.append(i)\n

\n

寻找要点：

\n

我们现在已经过滤了 3D 候选列表，现在是时候查找它们中是否有任何一个与我们的 2D 晶格匹配。为此，我们将使用三角形恒等式：如果两个三角形的所有边都相同，则它们是相同的。

\n

我们将如何做到这一点？

\n

我们已经有了 2D 和 3D 的距离矩阵。此外，我们还有起点及其可能的邻居。因此，对于每个这样的起点，我们将检查邻居的距离是否匹配。如果我们不进行任何优化，这将花费我们一些时间\xe2\x80\x94number_of_2d_points * number_or_3d_points。

\n

笔记：

\n

可能存在一些病态的边缘情况，这会发现太多匹配项。总有一种可能的优化：对距离进行排序，然后使用二分搜索进行匹配，以初始排序为代价，将每个点检查的时间从 n 减少到 log(n)。
\n令人烦恼的是，它会在索引周围打乱顺序，因此您需要记下它们并跟踪它们。

\n

算法：

\n

\n
• 采取起点
\n
• 通过距离检查找到其潜在邻居 (number_of_2d_points * number_of_3d_points)\n
  \n
  • 请注意，可能存在重复项。我们将尽快摆脱它们以防止组合爆炸
  \n
  • 一一检查（潜在的）邻居
  \n
  • 对于每个邻居\n
    \n
    • 检查与所有其他已检查邻居的距离
    \n
    • 如果全部匹配，则继续到下一个邻居
    \n
    • 否则丢弃该分支（如果有更多潜在邻居具有该距离，则尝试下一个，否则丢弃起点）
    \n
    \n
  \n
  • 如果我们匹配所有邻居，我们只是找到一组潜在的匹配点。（减去由于初始舍入而导致的不精确性，因此我们记下来）
  \n
  \n
\n

\n

笔记：

\n

我将尝试返回并提供代码并可能进行基准测试和优化。不过，我恐怕已经花了太多时间了，而且我还有一些职责要处理。

\n

现在怎么办？

\n

现在我们应该有一组（或多个）3D 匹配结构。如果我们没有找到，我们会降低精度并重试。如果我们找到多个，我们现在可以提高精度并查看（或计算）哪个候选者是最佳匹配\xe2\x80\x94这应该相当简单，因为我们已经在前一个点之间有了 1:1 映射步。（不过，跟踪指数会有点痛苦。）

\n

当我们获得最佳匹配时，最后一步要做的就是找到二维晶格方向。我们有：3D 点集（矩阵 P3d）、2D 晶格点坐标（矩阵 P2d）和晶格（矩阵 L2d）。
\n矩阵 P2d 和 P3d 的维数相同。
\n那么现在我们正在寻找这样一个矩阵 L3d，我们该如何实现呢？
\n我们需要找到一个矩阵 X，使得 X * P2d = P3d
\n然后我们做 X * L2d = L3d。

\n

笔记：

\n

有一种边缘情况，这将不起作用\xe2\x80\x94如果您找到一个结构，其中 2D 结构是 3D 结构的镜像，则不可能通过旋转来找到匹配项。如果你绘制它，应该很容易发现它。或者可能需要一些数学来找到这种情况。

\n

我们就完成了。

\n

如果我有更多时间，我会尝试用代码完成此工作并对其进行正确测试并可能进行优化。但不确定我什么时候能到达它。请告诉我这是否足够有帮助或者您是否需要我更深入地解释任何部分。
\n此外，如果有效，请告诉我进展如何:-)。

\n

编辑：

\n

回答您的问题：

\n

\n
1. 我已经对此有了一些想法并在上面列出了它们，但是现在由于尺寸差异我不确定。我将尝试简要总结它们。\n
   \n
   • 标准化两个晶格
   \n
   • 获取 3D 晶格并计算其边界顶点和边。
   \n
   • 取二维晶格的最长尺寸 2d_long（例如，连接彼此最远的两个顶点的线段的长度）
   \n
   • 计算以下边界框:\n
     \n
     • min_x = min_coord_x - 2d_long
     \n
     • max_x = max_coord_x + 2d_long
     \n
     • y 和 z 相同
     \n
     \n
   \n
   • 现在用晶格 3D 的副本填充此边界框。这应该是相当容易计算的。这应该足够了，希望不要太大。
   \n
   \n
\n
2. 好的，按细胞大小标准化，这应该不是问题，休息应该按照描述进行。它会影响精度，因此您可能只需要使用两位小数，随后会有更多匹配项
\n
3. 图意义上的邻居。\n
   \n
   • 我在 2D 晶格的顶点上创建一个完全完成的图（团），并指定边的权重等于长度。
   \n
   • 然后对于 3D 矩阵我做同样的事情，除了我丢弃不匹配的边\xe2\x80\x94e.g.，其长度与任何 2D 矩阵不同。
   \n
   • 最后我检查 3D 空间中是否存在匹配的子图。
   \n
   \n
\n
4. 您可以考虑角度匹配，但是如果您可以以足够的精度进行归一化，则不需要\xe2\x80\x94我们正在检查相同的三角形，并且这些三角形由于相同而具有匹配的角度。\n
   \n
   • 我们期望有多个匹配\xe2\x80\x94，事实上，在玩具示例中，有六个匹配的“起始”点。这并不重要\xe2\x80\x94我们将简单地检查它们并消除那些不会导致所需结构的内容。
   \n
   \n
\n

\n