非常大的图像中节省内存的线条拼接

13

背景

我处理来自合成孔径雷达卫星的非常大的数据集。这些可以看作是高动态范围灰度图像,每边约为10k个像素。

最近,我一直在开发单尺度变体的Lindeberg比例空间脊线检测算法的应用程序,用于在SAR图像中检测线性特征。这是对使用方向滤波器或使用霍夫变换的改进方法,因为它的计算成本比两者都要低。(我将在4月份的JURSE 2011上展示一些最新的结果,并且如果有帮助,我可以上传一份预印本)。

我目前使用的代码生成一个记录数组,每个像素对应一个记录,每个记录描述了与该像素右下方矩形内的一个脊线段相对应且由相邻像素限定的界限。

struct ridge_t { unsigned char top, left, bottom, right };
int rows, cols;
struct ridge_t *ridges;  /* An array of rows*cols ridge entries */

ridges中的一个条目包含一个山脊段,如果恰好有两个topleftrightbottom的值在0-128范围内。假设我有:

ridge_t entry;
entry.top = 25; entry.left = 255; entry.bottom = 255; entry.right = 76;

然后我可以找到山脊线段的起点 (x1,y1) 和终点 (x2,y2):
float x1, y1, x2, y2;
x1 = (float) col + (float) entry.top / 128.0;
y1 = (float) row;
x2 = (float) col + 1;
y2 = (float) row + (float) entry.right / 128.0;

当这些单独的脊段被渲染时,我得到了一个类似于这样的图像(一个更大图像的非常小的角落):

Rendered ridge segments

每个长曲线都由一系列小的脊线段渲染而成。
确定包含脊线段的相邻位置是否连接是微不足道的。如果我有位于(x, y)处的ridge1和位于(x+1, y)处的ridge2,则它们是同一条线的部分,如果0 <= ridge1.right <= 128且ridge2.left = ridge1.right
问题:理想情况下,我希望将所有的脊线段拼接成线,这样我就可以迭代在图像中找到的每一条线来应用进一步的计算。不幸的是,我发现很难找到一个既具有低复杂度又具有内存效率又适合多处理的算法(处理超大型图像时,这些都是重要考虑因素!)
我考虑过的一种方法是扫描整个图像,直到找到只有一个链接脊线段的脊,然后沿着结果的线走,标记线中的任何脊为已访问。但是,这对于多处理来说是不适用的,因为没有办法告诉是否有另一个线程从另一个方向(例如)走同一条线路,除非使用昂贵的锁定。
读者们有什么可能的方法建议吗?看起来像是过去有人已经想出了高效的方法...

1
对于这样一个写得如此出色的问题,给你点赞。 - user257111
4个回答

4

我不完全确定这是否正确,但我想让大家评论一下。首先,让我介绍一种无锁不交集算法,它将成为我提出的算法的重要组成部分。

无锁不交集算法

我假设您选择的CPU架构上存在两个指针大小的比较和交换操作。至少在x86和x64架构上可用。

该算法基本上与单线程情况下维基百科页面上描述的算法相同,但对于安全的无锁操作进行了一些修改。首先,我们需要将rank和parent元素都设置为指针大小,并在内存中对齐为2*sizeof(pointer),以便进行原子CAS操作。

Find()不需要更改;最坏的情况是在同时写入的情况下,路径压缩优化未能充分发挥作用。

Union()必须更改:

function Union(x, y)
  redo:
    x = Find(x)
    y = Find(y)
    if x == y
        return
    xSnap = AtomicRead(x) -- read both rank and pointer atomically
    ySnap = AtomicRead(y) -- this operation may be done using a CAS
    if (xSnap.parent != x || ySnap.parent != y)
        goto redo
    -- Ensure x has lower rank (meaning y will be the new root)
    if (xSnap.rank > ySnap.rank)
        swap(xSnap, ySnap)
        swap(x, y)
    -- if same rank, use pointer value as a fallback sort
    else if (xSnap.rank == ySnap.rank && x > y)
        swap(xSnap, ySnap)
        swap(x, y)
    yNew = ySnap
    yNew.rank = max(yNew.rank, xSnap.rank + 1)
    xNew = xSnap
    xNew.parent = y
    if (!CAS(y, ySnap, yNew))
      goto redo
    if (!CAS(x, xSnap, xNew))
      goto redo
    return

这应该是安全的,因为它永远不会形成循环,并且始终会得到正确的并集。我们可以通过观察以下事实来确认:
  • 首先,在终止之前,两个根之一将始终以指向另一个的父节点结束。因此,只要没有循环,合并就成功了。
  • 其次,排名始终增加。比较x和y的顺序后,我们知道在快照时x的排名低于y。为了形成循环,另一个线程需要先增加x的排名,然后合并x和y。但是,在写入x的父指针的CAS中,我们检查了排名是否发生了变化;因此,y的排名必须保持大于x。

在同时突变的情况下,y的排名可能会增加,然后由于冲突而返回重做。但是,这意味着y要么不再是根(在这种情况下,排名无关紧要),要么y的排名已经被另一个进程增加(在这种情况下,第二次尝试将不起作用,y将具有正确的排名)。

因此,不应该有形成循环的机会,这种无锁不相交集算法应该是安全的。

现在进入应用程序问题...

假设

我假设脊线段只能在其端点相交。如果不是这种情况,您将需要以某种方式更改第1阶段。

我还假设单个整数像素位置的共存足以连接脊线段。如果不是这样,您将需要更改第1阶段中的数组,以保存多个候选脊线段+不相交集对,并过滤出实际连接的那些。

此算法中使用的不相交集结构将在其结构中携带对线段的引用。在合并的情况下,我们任意选择两个记录的线段之一来表示集合。

第1阶段:本地线路识别

我们首先将地图分成扇区,每个扇区将作为单独的作业进行处理。可以在不同的线程中处理多个作业,但每个作业将仅由一个线程处理。如果脊线段跨越扇区,则将其分成两个部分,每个部分都属于一个扇区。

对于每个扇区,建立将像素位置映射到不相交集结构的数组。大部分此数组稍后将被丢弃,因此其内存要求不应太重。

我们首先处理扇区中的每个线段。我们首先选择一个不相交集合,表示线段所在的整条直线。我们首先在像素位置数组中查找每个端点,以查看是否已分配了不相交集合结构。如果其中一个端点已经在此数组中,则使用分配的不相交集合。如果两个端点都在数组中,则对不相交集合执行合并,并使用新的根作为我们的集合。否则,我们创建一个新的不相交集合,并将当前线段的引用与不相交集合结构相关联。然后,我们将我们的每个端点的新不相交集合的根写回到像素位置数组中。
这个过程对于扇区中的每个线段都要重复进行;最终,我们将通过不相交集合完全识别扇区内的所有线段。
请注意,由于不相交集合尚未在线程之间共享,因此还没有必要使用compare-and-swap操作;只需使用正常的单线程union-merge算法即可。由于在算法完成之前不释放任何不相交集合结构,因此可以从每个线程的bump分配器中进行分配,使内存分配(几乎)无锁且O(1)。
一旦一个扇区被完全处理,像素位置数组中的所有数据都将被丢弃;但是,与扇区边缘像素对应的数据将被复制到一个新数组中,并保留到下一阶段。
由于遍历整个图像是O(x*y),而不相交合并实际上是O(1),因此此操作为O(x*y),需要工作内存O(m+2*x*y/k+k^2)=O(x*y/k+k^2),其中t是扇区数,k是扇区的宽度,m是扇区中部分线段的数量(取决于线路跨越边界的频率,m可能会有很大变化,但永远不会超过线段数)。传递到下一个操作的内存为O(m+2*x*y/k)=O(x*y/k)
第2阶段:跨扇区合并
一旦所有扇区都已处理完毕,我们就转向合并跨越扇区的线路。对于每个扇区之间的边界,我们对跨越边界的线路执行无锁合并操作(即,在边界的每侧相邻像素已分配给线路集的地方)。
此操作的运行时间为O(x+y),消耗O(1)内存(但必须保留第1阶段的内存)。完成后,可以丢弃边缘数组。
第3阶段:收集线路
现在,我们在所有分配的不相交集合结构对象上执行多线程映射操作。我们首先跳过任何不是根的对象(即,当obj.parent != obj时)。然后,从代表线段开始,我们从那里向外移动,并收集和记录关于该线路的任何所需信息。我们可以确保每次只有一个线程查看任何给定的线路,因为相交的线路将最终进入同一不相交集合结构中。
这需要O(m)运行时间,并且内存使用量取决于您需要收集关于这些线段的哪些信息。

概要

总的来说,该算法的运行时间应为O(x*y),内存使用量为O(x*y/k + k^2)。调整k可以在第一阶段过程中在瞬时内存使用和传递到第二阶段的邻接数组和不相交集合结构之间进行权衡。

请注意,我实际上没有在真实世界中测试过这个算法的性能;也有可能我忽视了无锁不相交集合并算法中的并发问题。欢迎评论 :)


这看起来是一种非常有趣的方法;目前正在努力理解所有东西应该如何工作。 - Peter T.B. Brett
使用GCC的原子内建函数(http://gcc.gnu.org/onlinedocs/gcc-4.1.2/gcc/Atomic-Builtins.html)进行比较和交换操作是否合理? - Peter T.B. Brett
1
原来GCC将带有128位参数的__sync_bool_compare_and_swap()(顺便说一句,隐藏的__int128_t类型很有用)转换为调用__sync_bool_compare_and_swap_16(),而GCC没有提供实现!目前,我只在x86-64上运行此代码,因此我基于CMPXCHG16B提供了自己的实现。 - Peter T.B. Brett
1
我决定同时完成步骤1和2,然后将森林原地折叠成双向链表。结果非常惊人;在我的工作站上,对一个3500万像素的图像进行四进程线缝合只需要不到一秒钟的时间。感谢您的出色回答! - Peter T.B. Brett
嗨bdonlan,你能否发封电子邮件给我,我的邮箱是p.brett@surrey.ac.uk?谢谢。 - Peter T.B. Brett
显示剩余4条评论

2

你可以使用非泛化的Hough Transform算法。这种算法似乎在N x N网格数组上达到了令人印象深刻的O(N)时间复杂度(如果你有 ~10000x10000 SIMD数组,并且你的网格是 N x N - 注意:在你的情况下,N将是一条脊线结构,或者一个A x B脊线簇,而不是像素)。点击查看来源。更保守的(非核心)解决方案将复杂度列为O(kN^2),其中k = [-π/2, π]来源
然而,Hough变换确实需要一些相对较高的内存要求,空间复杂度将是O(kN),但如果您预先计算sin()和cos()并提供适当的查找表,则可以降低到O(k + N),这可能仍然太多,具体取决于您的N有多大...但我认为您无法再减少它了。
编辑:跨线程/内核/SIMD/进程线元素的问题是非常棘手的。我的第一反应是将网格细分为递归四叉树(取决于某个特定的公差),检查直接边缘并忽略所有边缘脊结构(实际上可以将这些标记为“潜在长线”并在整个分布式系统中共享);只对该特定四叉树内部的所有内容进行处理,并逐步向外移动。这是一个图形表示(绿色是第一次通过,红色是第二次通过等等)。然而,我的直觉告诉我,这是计算上昂贵的...

Passes


嗨,大卫!文献中已经有两种方法用于从SAR图像中进行线性特征提取: 方向滤波器组和霍夫变换。我正在进行的研究旨在尝试改进方向滤波器方法的精度,同时避免Hough变换中昂贵的优化步骤。令我不安的是,如果你从给定的脊结构开始,轻松地穿越整个图像,找到所有相连的脊段,但你无法确定另一个线程是否沿着相反的方向行走同一条路线。 - Peter T.B. Brett
@Peter:你的问题并不简单,我在上面的帖子中提供了一个暂时的解决方案,但我认为这可能是一种过于天真的方法。 - David Titarenco
嗨,David,谢谢你。那绝对是我考虑过的一种方法。我最担心的是需要单独跟踪工作单元内部行和工作单元间的行,但这绝对值得思考!我很想知道你对我刚刚添加为此问题提供的潜在方法的看法... - Peter T.B. Brett

0

好的,经过一番思考,我有一个建议,似乎太简单了,不够高效... 我希望能得到一些反馈,看看它是否合理!

1) 由于我可以轻松确定每个ridge_t山脊段是否连接到零个、一个或两个相邻的段,因此我可以适当地对每个山脊段进行着色(LINE_NONELINE_ENDLINE_MID)。这可以很容易地并行完成,因为不存在竞争条件。

2) 着色完成后:

for each `LINE_END` ridge segment X found:
    traverse line until another `LINE_END` ridge segment Y found
    if X is earlier in memory than Y:
        change X to `LINE_START`
    else:
        change Y to `LINE_START`

这也是没有竞态条件的,因为即使两个线程同时遍历同一行,它们也会做出相同的更改。

3)现在图像中的每一行都将有一个端点被标记为LINE_START。可以在单个线程中定位并打包这些线条到一个更方便的结构中,而无需进行任何查找以查看该线条是否已被访问。

我应该考虑在第2步中收集诸如线条长度之类的统计信息,以帮助最终重新打包...

我有没有忽略任何陷阱?

编辑:明显的问题是我最终需要走两次线路,一次是为了定位RIDGE_START,一次是为了进行最终的重新打包,导致计算效率低下。尽管如此,从存储和计算时间的角度来看,它仍然似乎是O(N),这是一个好迹象...


如果直线相交会发生什么? - bdonlan
啊,糟糕。好发现。我想还需要循环检测。 - Peter T.B. Brett

0
如果山脊足够清晰,断点只有几个像素,那么查找线条/OCR所需的标准膨胀-查找邻居-侵蚀步骤应该可以解决。
从许多段连接更长的轮廓,并知道何时创建颈部或何时创建单独的岛屿要复杂得多。

因为我知道脊的位置具有亚像素精度,并且因为我可以轻松测试两个相邻的脊段是否连接,所以像膨胀/邻居/侵蚀这样的有损算法并不是我要寻找的。事实上,使用脊检测方法的整个原因是为了避免传统线查找/OCR方法的缺点。是的,我知道我所要求的要复杂得多。;-) - Peter T.B. Brett
好的 - 我唯一的经验是从切割TIN创建轮廓,其中你会得到许多不连续的直线段,但顺序混乱。 - Martin Beckett

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接