具有周期性边界条件的 3D 连接组件

Question

我正在尝试识别地球上的连接特征（即球体上的时空）。 cc3d 包已经完成了 90% 的工作，但我正在努力处理日期边界（即三个维度之一的周期性边界条件）。

在 2D 地图上，我的方法的问题变得明显（注意南极附近 0 经度周围的连接购买不同标记的区域）：

出现这种情况是因为数据是在 0 到 360 的经度上定义的（虽然我在这里显示了从 -180 到 180 的经度，以使问题更加明显）。

我针对此类问题的两种典型解决方案不起作用：

• 将反经络翻转到太平洋只会转移问题，因此没有帮助。
• 在右侧边界连接数据副本也无济于事，因为它会导致原始左侧和右侧粘贴数据之间的标签不明确

MWE

问题在二维上的分解应该如下：

import cc3d
import numpy as np
import matplotlib.pyplot as plt

arr = np.sin(np.random.RandomState(0).randn(100).reshape(10, 10)) > .3
labels = cc3d.connected_components(arr)
map = plt.pcolormesh(labels, vmin=.5)
map.cmap.set_under('none')

此处右上角的黄色结构应连接到顶部结构的其余部分，底部的两个结构也是如此。请记住，任何有用的解决方案也应该适用于 3 维连接特征。

任何有关如何解决此问题的帮助都将不胜感激。

Answer 1

2024 年 4 月更新： 周期性边界条件选项现已直接在 cc3d 中实现（请参阅下面 SapphireSun 的评论）。

---

好吧，我已经修复了输出。它非常慢，需要运行

connected_components

两次，但对我来说很有效。我把它放在这里，以防以后对任何人有帮助。

我将在这里修复 MWE，但它对于完整的 3D 时空数据同样有效（每天独立）。 我正在使用的数据的约定是经度是最后一个维度，因此（时间、纬度、经度），但它绘制在 x 轴上，因此我对所有函数调用执行

swapaxes

MWE 示例。

功能

def plot(arr):
    map = plt.pcolormesh(arr, vmin=.5)
    map.cmap.set_under('none')    
 
    for xx in range(len(arr)):
        for yy in range(len(arr)):
            plt.text(yy+.5, xx+.5, arr[xx, yy],
                     ha="center", va="center", color="k")


def cut_stitch(arr):  # equivalent to switching the antimeridian [0, 360[ -> [-180, 180[
    idx = len(arr) // 2
    return np.concatenate([arr[idx:], arr[:idx]], axis=0)


def fix_date_border(arr, arr_r):
    mask = arr == 0

    mapping = {}
    for key, value in zip(arr_r[~mask], arr[~mask]):
        try:
            mapping[key].add(value)
        except KeyError:
            mapping[key] = {value}  
    # if there is only one unique value per key its a one-to-one mapping and no problem
    conflicts = [list(values) for values in mapping.values() if len(values) != 1]
    
    print(mapping)  # DEBUG
    
    for conflict in conflicts:
        idx = np.any([arr == cc for cc in conflict[1:]], axis=0)
        arr[idx] = conflict[0]  # reassign in original (i.e., first) array
        
    return arr

在重新排列的字段上第二次运行组件标签，然后再次重新排列：

arr2 = cut_stitch(arr.swapaxes(0, 1)).swapaxes(0, 1)
labels2 = cut_stitch(cc3d.connected_components(arr2).swapaxes(0, 1)).swapaxes(0, 1)
plot(labels)  # original
plot(labels2)

修复标签：

labels_fixed = fix_date_border(labels.swapaxes(0, 1), labels2.swapaxes(0, 1)).swapaxes(0, 1)
# {2: {1, 2}, 7: {8, 7}, 5: {4}, 3: {3}, 1: {1}, 4: {5}, 8: {7}, 6: {6}}
plot(labels_fixed, labels=True)

请注意，固定数组中的标签将不再严格连续！