骨架化密集标记的图像体积。

kimimaro的Python项目详细描述


构建状态pypi版本

Kimimaro:骨架化密集标记图像

使用teasar导出的方法快速骨架化二维和三维numpy阵列中的所有非零标签。返回的骨架列表的格式为云卷

在3.7GHz英特尔i7处理器上,此软件包在一分钟内处理了一个512x512x100的卷,包含333个标签。它在8到13分钟内处理了一个512x512x512的卷,其中包含2124个标签(取决于是否设置了fix\u branching

用kimimaro标记的密集体积骨架

pip安装

< EM >需要C++编译器。</P>

sudo apt-get install python3-dev g++
pip3 install numpy
pip3 install kimimaro 

将来,我们可能会创建一个完全二进制分布。

示例

用kimimaro标记的密集体积骨架
图2:512x512x512密集标记卷上的内存使用情况

图2显示了将Kimimaro 0.5.2应用于包含2124个连接组件的Connectomics数据集的512x512x512断开标签时所需的内存使用和处理时间(约390秒,约6.5分钟)。描述了算法的不同部分。严重的是,序言运行大约半分钟,骷髅化大约六分钟,并在几秒钟内完成。峰值内存使用量约为4.1 GB。下面的代码用于处理标签。由于fix_borders和max_path的组合,glia的处理在中被截断。

Kimimaro已经走了很长一段路。版本0.2.1花费了15分钟,在同一数据集上的前导码运行时间是前者的两倍。

# LISTING 1: Producing Skeletons from a labeled image.importkimimarolabels=np.load(...)skels=kimimaro.skeletonize(labels,teasar_params={'scale':4,'const':500,# physical units'pdrf_exponent':4,'pdrf_scale':100000,'soma_detection_threshold':1100,# physical units'soma_acceptance_threshold':3500,# physical units'soma_invalidation_scale':1.0,'soma_invalidation_const':300,# physical units'max_paths':50,# default None},# object_ids=[ ... ], # process only the specified labelsdust_threshold=1000,# skip connected components with fewer than this many voxelsanisotropy=(16,16,40),# default Truefix_branching=True,# default Truefix_borders=True,# default Trueprogress=True,# default False, show progress barparallel=1,# <= 0 all cpu, 1 single process, 2+ multiprocess)# LISTING 2: Combining skeletons produced from #            adjacent or overlapping images.importkimimarofromcloudvolumeimportPrecomputedSkeletonskels=...# a set of skeletons produced from the same label idskel=PrecomputedSkeleton.simple_merge(skels).consolidate()skel=kimimaro.postprocess(skel,dust_threshold=1000,# physical unitstick_threshold=3500# physical units)

调整kimimaro.skeletonie参数

该算法通过在三维物体上寻找一个根点,然后通过dijksta的最短路径算法通过一个惩罚域将路径串行跟踪到最远的不可见点。每次经过后,都会有一个球体(实际上是一个外切立方体)围绕当前路径中的每个顶点展开,该球体会标记所访问对象的一部分。

有关骨骼化过程基本知识的可视化教程,请参阅以下wiki文章:一篇关于teasar骨骼化的图片指南

更详细的信息,阅读下面的内容茶叶纸(尽管我们德维亚在一些地方)。〔1〕

比例常量

通常,需要调整的最重要的参数是scaleconst这两个参数根据以下方程式控制这个失效球体的半径:r(x,y,z)=scale*dbf(x,y,z)+const其中尺寸是物理的(例如纳米,即校正各向异性)。dbf(x,y,z)是从该点到形状边界的物理距离。

查看这篇wiki文章来帮助你完善你的直觉。 各向异性

表示每个体素的物理维度。例如,可以用电子显微镜以每像素4nm×4nm的速度扫描连接组学数据集,并将其堆叠成40nm厚的切片。即:<编码>各向异性=(4,4,40)。只要一致,就可以使用任何单位。

灰尘阈值

此阈值剔除小于此多个体素的连接组件。

最大路径

限制可为给定标签绘制的路径数。某些对当前分析可能不重要的细胞,如胶质细胞,处理起来可能很昂贵,可能会提前中止。

pdrf_标度pdrf_指数

pdrf_scalepdrf_exponent表示惩罚方程的参数,该惩罚方程采用欧几里德距离场(d)并对其进行扩展,以便更接近边界的切割是非常不利的,从而使dijkstra走的路径更居中。

pr=pdrf_标度*(1-d/max(dpdrf_指数+(方向梯度<;1.0)。

默认设置应该可以很好地工作,但是在大的各向异性或海绵状形态下,可能需要对其进行调整。如果你看到骨架在一个大的区域内乱糟糟的,那可能是浮点精度的崩溃。

soma_acceptance_thresholdsoma_detection_threshold

我们特别处理somas,因为它们没有管状几何体,而是应该以中心和辐条的方式表示。soma_acceptance_threshold是物理半径(例如,以纳米为单位),在此半径之外,我们将图像的连接组件分类为包含soma。距离变换的输出被标签上的空洞所抑制,这些空洞通常是由somata上的分割算法产生的。我们可以填充它们,但是我们使用的孔填充算法很慢,所以我们只希望偶尔应用它。因此,我们设置了一个较低的阈值,即soma_acceptance_threshold,超过该阈值,我们将填充孔并重新测试soma。

soma_invalization_scalesoma_invalization_const

一旦我们将一个区域分类为soma,我们就将骨架化算法的根固定在离边界最大距离的一个点上(通常只有一个)。然后,我们将该点周围的所有体素标记为已访问,其球形半径由r(x,y,z)=soma_invalization_scale*dbf(x,y,z)+soma_invalization_const描述,其中dbf(x,y,z)是距该点形状边界的物理距离。如果做得正确,这可以防止骨骼被绘制到胞体的边界,而是将骨骼主要拉到从细胞体延伸的过程中。

修复边框

此功能使连接不适合RAM的相邻图像卷的骨架变得更容易。如果启用,骨骼将被确定地绘制到形状与边界接触的每个位置的二维接触区域的近似中心。这会影响操作的性能取决于触点的形状和数量,可以是正触点,也可以是负触点。

修复分支

你可能永远不会想禁用它,但基本茶是臭名昭著的叉在分支点过早的骨架。此选项使您可以优先选择在更合理的位置叉取,但会导致严重的性能损失。

进度

骨架化阶段开始后显示进度条。

并行

使用一个处理器池来更快地实现骨架化。每个进程可分配的任务都是一个连接组件的骨架化(因此,对于一个需要很长时间才能骨架化的标签没有帮助)。此选项还影响初始欧氏距离变换的速度,该变换是并行启用的,是前导码中最昂贵的部分(如下所述)。

性能提示

动机

连接组学领域通常产生大量密集标记的神经组织。骨骼是二维或三维物体的一维表示。它们有许多用途,其中一些是神经元的可视化、计算全局拓扑特征、快速测量对象之间的电距离,以及将树结构强加于神经元(对计算和用户界面有用)。有几种方法可以计算骨骼,也有几种方法可以定义它们[4]。经过一些实验,我们发现teasar[1]方法给出了相当好的结果。其他方法包括拓扑细化("洋葱剥皮")和寻找由最大内接球体描述的中心线。Seung Lab的校友Ignacio Arganda Carreras为斐济编写了一个拓扑细化插件,名为"Skeletonize3D"。

在connectomics领域中有几种teasar的实现方法[3][5],但是人们通常都知道teasar的实现是缓慢的,并且可以使用几十亿字节的内存。我们的目标是在一个花瓣体素比例的图像中快速骨架化所有标签,这清楚地表明现有的稀疏实现是不切实际的。在将稀疏方法应用于云管道时,我们注意到在重复计算欧氏距离变换时存在效率低下的问题。rm(edt),连接分量算法的重复计算,在dijkstra算法使用的图的构造中,边由体素之间的空间关系表示,在存储成本中,体素数量是二次方的,用于表示图像中隐含的图形、用于表示相对较小的剪切块的不必要的大数据类型以及重叠区域的重复下载。我们还发现,teasar的"滚动失效球"的简单实现不必要地以骨架路径长度的二次方来重新评估大量体素。

我们进一步发现,edt的商品实现只支持二进制映像。我们无法找到任何可用的Python或C++库来执行图像上的Dijkstra最短路径。图像的连接组件算法的商品实现仅支持二进制图像。因此,设计了几个库来弥补这些缺陷(参见相关项目)。

为什么是提萨?

teasar:m.sato等人在2000年发表的一篇论文[1]中提出了一种用于精确和健壮骨骼的树结构提取算法,该算法是将二维和三维结构转换为嵌入该高维中的一维"骨架"的一系列算法的一员。人们可以将骨骼简化为从二值图像中提取棒状图形。这个问题比看起来更难。在画这样一幅画时,必须考虑不同的情况。例如,香蕉的棒状图可能只是一个弯曲的中心线,而甜甜圈的图可能是一个闭环。在我们分析神经元的案例中,有时我们希望骨骼包括棘突,通常与突触相连的树突的短突起,有时我们只需要描述一个轴突主干的长度。

此外,数据质量问题也可能具有挑战性。如果一个人正在对一个甜甜圈的二维图像进行骨架化,但是角度与圆环的正交轴有足够的偏差,那么甚至有可能精确地执行这项任务吗?在三维情况下,如果一个神经元的标记出现断裂或合并,算法是否能正常工作?这些问题在手动和自动图像精液中都很常见。

在我们的各向异性体素标记的骨架化神经元问题领域中,我们选择的算法应该产生树结构,根据具体情况处理精细或粗略的细节提取,处理体素各向异性,并且在cpu和内存使用方面具有合理的效率。茶满足了这些标准。值得注意的是,teasar并不能保证骨架在形状中的中心位置,但它做出了努力。基本的teasar算法是已知的在转弯处抄近路和过早分支。2001年,teasar团队的成员在第204页第4.2.2节中描述了一种减少早期分支问题的方法。〔2〕

teasar衍生算法

我们实现了teasar算法,但是为了提高路径中心性,提高性能,处理膨胀的细胞体,并实现对大图像的有效分块评估,我们与现有算法有了一些偏差。我们选择不实现[2]中的梯度向量场步骤,因为我们的实现已经相当快了。这篇论文声称输入体素减少了70-85%,因此值得研究。

为了处理包含许多标签的图像,我们的一般策略是尽可能多地执行操作,以便在一次传递中处理所有标签。我们实现的一些组件算法(例如,连接组件、欧几里德距离变换)每次通过可能需要几秒钟,因此,重要的是它们不能运行数百或数千次。本PA的大部分工程贡献ckage取决于这些操作的效率,这些操作将运行时间从小时减少到分钟。

给定一个三维标记体素阵列,i,n>;=0个标记,有序的三重描述体素各向异性a,我们的算法可以分为三个阶段,按照这个顺序进行预处理、骨架化和定形。

一、序言

前导码获取包含n标签的3d图像,并有效地生成骨骼化阶段所需的连接组件、距离变换和边界框。

  1. 为了提高性能,如果n为0,则返回一组空骨架。
  2. 标记m个连接的组件,icc,ofi
  3. 要节省内存,请按从1到m的顺序对连接的组件重新编号。将新图像的数据类型调整为包含m并覆盖icc的最小uint类型。
  4. 生成重新编号的icci的映射,以便以后为骨骼分配有意义的标签,并删除i以节省内存。
  5. 计算e,给定aicc的多标签各向异性欧氏距离变换。e将所有隔行扫描边视为变换边,而不是图像的边界。黑色像素视为背景。
  6. icc中收集唯一标签的列表,lcc和基于它们所代表的体素数量来处理哪些标签以去除"灰尘"。
  7. 在一次过程中,计算与lcc中每个标签对应的边界框列表b。 < > >

    Ⅱ。骨骼化

    在此阶段,我们从每个连接的组件标签中提取树结构的骨架。下面,我们引用序言中定义的变量。为了清楚起见,我们省略了soma特定的处理,并保持fix\u branching=true

    对于lccb中的每个标签l

    1. 提取il,使用blicc中紧密包围l的裁剪二值图像
    2. 使用ilble中提取elell的裁剪紧密封闭的edt。这比为每个二进制图像重新计算EDT快得多。
    3. 找到任意前景体素并使用该点作为源,计算il的各向异性欧几里德距离场。最大值的坐标现在是"根"r
    4. r计算欧氏距离场,并将其保存为与根场的距离dr
    5. 计算从根字段pr=pdrf评分*((1-el/max(el)^pdrf指数)+dr/max(dr)。
    6. il包含前景体素:
      1. 将目标坐标t确定为dr中距离最大的前景体素。
      2. 考虑到pr中的体素值作为边权重,绘制从rt的最短路径。
      3. 对于p中的每个顶点v,扩展物理边长为scale*elv)+const的无效立方体,并将与这些立方体重叠的il中的任何前景像素转换为背景像素。
      4. (仅当p中的每个顶点坐标v固定分支=真)时,设置prv)=0。
      5. p附加到此标签的路径列表。
      6. < > >
      7. 使用el,提取到骨架中每个顶点表示的最近边界的距离。
      8. 对于从il中提取的每个原始骨架,将顶点平移bl,以更正由裁剪操作引起的平移。
      9. 将顶点乘以各向异性a将它们放置在物理空间中。
      10. < > >

        如果考虑soma处理,我们修改根(r)搜索过程如下:

        1. 如果max(el)>;soma_detection_threshold..
        2. il中填充拓扑孔。SOMA是一个大的区域,通常有来自不完善的自动标记方法的灰尘。
        3. 从这个清理过的图像中重新计算el
        4. 如果max(el)>;soma_acceptance_threshold,则转到soma处理模式。
        5. 如果处于SOMA处理模式,则继续,否则转到上面算法中的步骤3。
        6. r设置为与max对应的坐标(el
        7. 创建一个物理半径为soma_invalization嫒u scale*max(el)+soma_invalization嫒u const的无效球体,并从其中包含的il中擦除前景体素。这有助于防止在SOMA上绘制出错误的路径。
        8. 从上述算法的步骤4继续。
        9. < > >

          三、定稿

          在最后阶段,我们将不同连接的组件骨架聚集成单个骨架,并将它们的标签分配给输入图像。这一步是人为的,与它的实现与骨架化的混合方式相比,它是分开的,但在概念上是分开的。

          与teasar的偏差

          有几个地方我们采取了不同于teasar作者所要求的方法。

          将DAF用于目标,PDRF用于寻路

          原始teasar算法定义了与根体素场(pdrf,pr以上)的惩罚距离为:

          PDRF = 5000 * (1 - DBF / max(DBF))^16 + DAF
          

          dbf是到边界场的距离(el以上),daf是到任何体素场的距离(dr以上)。

          我们发现,daf的加入往往会干扰从中心线开始的骨架路径,而这更好地由反向dbf单独描述。我们还发现修改常数和指数有助于调整转弯行为。最初,我们完全去掉了pdrf中daf的添加,但这引入了另一种问题。pdrf的指数运算导致浮点值在很大的开放空间中崩溃。这使得骨骼在追踪由浮点错误描述的路径时变得疯狂。

          daf在根体素和目标体素之间提供了一个非常有用的渐变,我们只是不希望这个渐变将路径偏离中心线。因此,鉴于pdrf基场非常大,我们加入了归一化daf,它刚好足以克服浮点误差,并在宽管和凸起处提供方向。

          原始文件还要求使用最大(pdrf)前景值选择目标。然而,这有点奇怪,因为pdrf值由边界效应而不是纯距离度量支配。因此,我们从最大值(daf)中选择目标。

          零加权先前路径(fix_branching=true

          2001年的骨架化论文[2]呼吁通过使用已经计算出的路径顶点作为场源计算daf来校正早期分叉。这使得dijkstra的算法可以免费追踪现有的路径,并在离目标较近的地方偏离它。

          由于我们已经强烈地不强调daf在dijkstra路径发现中的作用,因此不需要计算这个字段,我们只需要沿着e的路径将pdrf设置为零。为达到这个效果而存在的骨骼。这为我们节省了每个路径的昂贵的重复daf计算。

          然而,由于我们一直在计算一个dijkstra"父域",它记录了从每个前景体素到根的最短路径,所以采用这种方法仍然要付出很大的代价。然后我们使用这个保存的结果快速计算所有路径。然而,由于这种零加权修改使得后续计算依赖于先前的计算,我们需要为每条路径重新计算dijkstra算法。

          非重叠分块处理(fix_borders=true

          在处理大型卷时,批量生成骨架的一个明智方法是对卷进行分块,独立处理分块,并在最后合并生成的骨架片段。然而,这是复杂的"边缘效应"导致的上下文丢失,这使得它不可能期望由相邻块产生的骨架片段的端点对齐。相比之下,很容易连接网格碎片,因为网格碎片边缘的顶点位于给定一个重叠像素的可预测的相同位置。

          以前,我们使用50%的重叠来连接相邻的骨架片段,这将大体积骨架化的计算成本增加了8倍。但是,如果我们可以强制骨骼位于边界上的可预测位置,我们可以使用单像素重叠并复制简单的网格连接方法。作为一种(不正确但有用的)直觉,考虑计算每个边界平面上每个连接组件的质心,并将其添加为所需的路径目标。这将保证平面的两边连接在同一个像素上。但是,质心可能不在非凸面外壳的内部,因此我们必须更加复杂,并在形状内部选择一些真实的点。

          为此,我们再次重新调整欧氏距离变换的用途,并将其应用于连接组件的六个平面中的每一个,并选择最大值作为强制目标。这对于许多类型的对象都很有效,这些对象与单个平面接触并且具有单个最大值。但是,我们必须处理具有多个最大值的长方体和形状的角。

          为了处理与长方体多个侧面接触的形状,我们只需将目标分配给所有连接的组件。如果这在后处理中引入了一个循环,那么我们已经有了在火成岩中处理它的循环移除代码。如果它引入了一些无用的小附件,我们也有代码来处理这个问题。

          如果形状具有多个距离变换最大值,则在不需要在可能在不同机器上不同时间运行的空间相邻任务之间通信的情况下选择相同的像素是很重要的。此外,相邻任务上的同一平面会翻转坐标系。一种简单的方法可能是在其中一个坐标帧中选择具有最小x和y(或其他基于坐标的标准)的坐标,但这需要跟踪所有六个平面上的翻转,而且很烦人。取而代之的是,我们使用了一系列基于无坐标拓扑的过滤器,这些过滤器既有趣,又省力,而且选择了一些外观合理的东西。对这种方法的一个有效的批评是,它会在一个完全对称的物体上失败,但这些物体在生物数据中是罕见的。

          我们应用一系列过滤器,并根据通过的第一个过滤器选择点:

          1. 最接近当前标签质心的体素。
          2. 最接近图像平面质心的体素。
          3. 最接近飞机的一角。
          4. 最接近平面边缘。
          5. 先前找到的最大值。
          6. < > >

            重要的是,filter 1必须基于标签的形状,以便使凸面外壳的扭结最小化。例如,最初我们只使用2到5个过滤器,但这导致远离大块中心的神经突骨骼突然在区块边界处向区块中心移动。

            滚动失效立方体

            最初的teasar论文呼吁在步骤6(iii)中使用"滚动无效球"来擦除前景体素。这个球的一个简单实现非常昂贵,因为路径中的每个体素都需要自己的球,而且这些体素中有许多重叠。在某些情况下,可能需要对从根到目标路径上的每个体素无意义地重新评估整个体积。虽然在最坏的情况下可能会出现特殊情况,但在更常见的一般情况下,会花费大量的重复工作。

            因此,我们应用了一种使用拓扑线索的算法来执行线性时间内的失效操作。为了实现的简单,我们用立方体代替了球体。函数名roll\u invalization\u cube旨在引起这种尴尬,尽管它似乎并不重要。

            两遍算法如下。给定一个二值图像i、一个骨架s和一组顶点v

            1. bv为一组边界框,这些边界框将茶叶纸所示的球体标记出来。
            2. 分配一个3d有符号整数数组,t,表示拓扑的i的大小和维数。t最初设置为全零。
            3. 对于每个bv
            4. 对于沿x轴的bv左边界上的所有点p设置t(p)+=1。
            5. 将t(p)-=1设置为沿x轴的所有点的右边界上的ponbv
            6. 计算边界框bglobal表示所有bv的并集
            7. 从yz平面开始,每行b全局的点p沿x轴移动。
            8. 设置整数着色=0
            9. 在每个索引处,着色+=t(p)
            10. 如果着色>;0或t(p)为非零(我们在左边缘),则我们位于无效立方体中,并开始将前景体素转换为背景体素。
            11. < > >

              相关项目

              为了使这个模块成为可能,必须对几个经典算法进行特殊调整。

              1. edt:支持euclidean距离变换实现的单通道、多标签各向异性。
              2. dijkstra3d:dijkstra在26个连接的3d图像上定义的最短路径算法。这避免了边生成的时间开销和图形表示的内存浪费。
              3. connected-components-3d:在26个具有多个标签的连接的3d图像上定义的连接组件实现。
              4. fastremap:允许从3d数组中的1高速重新编号标签,以减少不必要的大32位和64位标签造成的内存消耗。
              5. < > >

                此模块最初设计用于CloudVolume和火成岩。

                1. cloudvolume:用于读取和写入神经组织、网格和骨骼的petascale分块图像的无服务器客户端。
                2. 火成岩:可视化连接组学数据集的分布式计算。
                3. < > >

                  此软件包中使用的一些teasar修改首先由alex bae演示。

                  1. 骨架化:针对稀疏标签的修改teasar的python实现。
                  2. < > >

                    学分

                    亚历克斯·贝克·戴夫介绍了前体骨架化包和我们在该包中使用的teasar的一些改进。亚历克斯还开发了后处理方法,用于使用50%重叠拼接骨骼。will silversmith将这些技术用于大规模生产,改进了一些基本的算法来同时处理数千个标签,并将它们重写到kimimaro包中。将添加涓流daf,零加权先前探索的路径,并固定边界的算法。forrest collman增加了参数灵活性,并帮助调整daf计算性能。Sven Dorkenwald和Forrest都提供了有益的讨论和反馈。

                    参考文献

                    1. 佐藤先生、苦味先生、本德先生、考夫曼先生和中岛先生。"teasar:用于精确和健壮骨骼的树结构提取算法。PROC第八届太平洋计算机图形学与应用会议。10月2000日。doi:10.1109/pccga.2000.883951(链接
                    2. I.苦,A.E.考夫曼和M.佐藤。"惩罚距离体积骨架算法"。ieee可视化与计算机图形学汇刊第7卷,国际空间站。2001年7月至9月3日。doi:10.1109/2945.942688(链接
                    3. 赵先生,S广场。"通过果蝇延髓内轴突定位自动识别神经元类型。9月2014日。arxiv:1409.1892[q-bio.nc](链接
                    4. A.Tagliasacchi、T.Delame、M.Spagnuolo、N.Amenta、A.Telea。"三维骨骼:最新的报告"。2016年5月。计算机图形学论坛。第35卷,国际空间站。2。doi:10.1111/cgf.12865(链接
                    5. P.Li、L.Lindsey、M.Januszewski、Z.Zheng、A.Bates、I.Taisz、M.Tyka、M.Nichols、F.Li、E.Perlman、J.Maitin Shepard、T.Blakely、L.Leavitt、G.Jefferis、D.Bock、V.Jain。"利用洪水充盈网络和局部重新调整,对果蝇脑的连续切片进行自动重建"。2019年4月。比奥希夫doi:10.1101/605634(链接
                    6. < > >

                      欢迎加入QQ群-->: 979659372 Python中文网_新手群

                      推荐PyPI第三方库


热门话题
在java中构建rest api时,使用servlet有哪些优点   java如何在Get上初始化@DateTimeFormat参数?   安卓 Java图像阵列   java如何在使用Hibernate访问器时忽略pmd LocalVariableCouldBeFinal   java安卓 If语句失败   java在Joda时间内设置默认年份   java Clojure安装:启动Clojure失败   com包之间的java冲突。fasterxml。杰克逊。数据绑定。jsonschema和com。fasterxml。杰克逊。数据绑定。jsonSchema   java数组返回null。。生成随机字   java对main方法的有趣调用   java将实体从onetomany关系重构为onetoone关系(它有自己的onetomany,便于更新)有意义吗?   json序列化/反序列化期间的java句柄循环引用   java在没有应用服务器的情况下,如何伪造数据源的jndi查找   java JDBC将转到错误的地址   从Java执行sonar Embedded Runner后,内存不足无法删除persistit_tempvol文件   java如何以Json格式显示包含员工和多个地址的地图   java在Android中同时从多个Firebase数据库引用中检索数据