骨架化密集标记的图像体积。
kimimaro的Python项目详细描述
Kimimaro:骨架化密集标记图像
使用teasar导出的方法快速骨架化二维和三维numpy阵列中的所有非零标签。返回的骨架列表的格式为云卷
在3.7GHz英特尔i7处理器上,此软件包在一分钟内处理了一个512x512x100的卷,包含333个标签。它在8到13分钟内处理了一个512x512x512的卷,其中包含2124个标签(取决于是否设置了fix\u branching
。
pip
安装
< EM >需要C++编译器。</P>
sudo apt-get install python3-dev g++ pip3 install numpy pip3 install 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〕
比例
和常量
通常,需要调整的最重要的参数是scale
和const
这两个参数根据以下方程式控制这个失效球体的半径:r(x,y,z)=scale*dbf(x,y,z)+const其中尺寸是物理的(例如纳米,即校正各向异性)。dbf(x,y,z)
是从该点到形状边界的物理距离。
查看这篇wiki文章来帮助你完善你的直觉。
各向异性
表示每个体素的物理维度。例如,可以用电子显微镜以每像素4nm×4nm的速度扫描连接组学数据集,并将其堆叠成40nm厚的切片。即:<编码>各向异性=(4,4,40)。只要一致,就可以使用任何单位。 此阈值剔除小于此多个体素的连接组件。 限制可为给定标签绘制的路径数。某些对当前分析可能不重要的细胞,如胶质细胞,处理起来可能很昂贵,可能会提前中止。 pr= 默认设置应该可以很好地工作,但是在大的各向异性或海绵状形态下,可能需要对其进行调整。如果你看到骨架在一个大的区域内乱糟糟的,那可能是浮点精度的崩溃。 我们特别处理somas,因为它们没有管状几何体,而是应该以中心和辐条的方式表示。 一旦我们将一个区域分类为soma,我们就将骨架化算法的根固定在离边界最大距离的一个点上(通常只有一个)。然后,我们将该点周围的所有体素标记为已访问,其球形半径由 此功能使连接不适合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算法,但是为了提高路径中心性,提高性能,处理膨胀的细胞体,并实现对大图像的有效分块评估,我们与现有算法有了一些偏差。我们选择不实现[2]中的梯度向量场步骤,因为我们的实现已经相当快了。这篇论文声称输入体素减少了70-85%,因此值得研究。 为了处理包含许多标签的图像,我们的一般策略是尽可能多地执行操作,以便在一次传递中处理所有标签。我们实现的一些组件算法(例如,连接组件、欧几里德距离变换)每次通过可能需要几秒钟,因此,重要的是它们不能运行数百或数千次。本PA的大部分工程贡献ckage取决于这些操作的效率,这些操作将运行时间从小时减少到分钟。 给定一个三维标记体素阵列,i,n>;=0个标记,有序的三重描述体素各向异性a,我们的算法可以分为三个阶段,按照这个顺序进行预处理、骨架化和定形。 前导码获取包含n标签的3d图像,并有效地生成骨骼化阶段所需的连接组件、距离变换和边界框。 在此阶段,我们从每个连接的组件标签中提取树结构的骨架。下面,我们引用序言中定义的变量。为了清楚起见,我们省略了soma特定的处理,并保持 对于lcc和b中的每个标签l… 如果考虑soma处理,我们修改根(r)搜索过程如下: 在最后阶段,我们将不同连接的组件骨架聚集成单个骨架,并将它们的标签分配给输入图像。这一步是人为的,与它的实现与骨架化的混合方式相比,它是分开的,但在概念上是分开的。 有几个地方我们采取了不同于teasar作者所要求的方法。 原始teasar算法定义了与根体素场(pdrf,pr以上)的惩罚距离为: dbf是到边界场的距离(el以上),daf是到任何体素场的距离(dr以上)。 我们发现,daf的加入往往会干扰从中心线开始的骨架路径,而这更好地由反向dbf单独描述。我们还发现修改常数和指数有助于调整转弯行为。最初,我们完全去掉了pdrf中daf的添加,但这引入了另一种问题。pdrf的指数运算导致浮点值在很大的开放空间中崩溃。这使得骨骼在追踪由浮点错误描述的路径时变得疯狂。 daf在根体素和目标体素之间提供了一个非常有用的渐变,我们只是不希望这个渐变将路径偏离中心线。因此,鉴于pdrf基场非常大,我们加入了归一化daf,它刚好足以克服浮点误差,并在宽管和凸起处提供方向。 原始文件还要求使用最大(pdrf)前景值选择目标。然而,这有点奇怪,因为pdrf值由边界效应而不是纯距离度量支配。因此,我们从最大值(daf)中选择目标。 2001年的骨架化论文[2]呼吁通过使用已经计算出的路径顶点作为场源计算daf来校正早期分叉。这使得dijkstra的算法可以免费追踪现有的路径,并在离目标较近的地方偏离它。 由于我们已经强烈地不强调daf在dijkstra路径发现中的作用,因此不需要计算这个字段,我们只需要沿着e的路径将pdrf设置为零。为达到这个效果而存在的骨骼。这为我们节省了每个路径的昂贵的重复daf计算。 然而,由于我们一直在计算一个dijkstra"父域",它记录了从每个前景体素到根的最短路径,所以采用这种方法仍然要付出很大的代价。然后我们使用这个保存的结果快速计算所有路径。然而,由于这种零加权修改使得后续计算依赖于先前的计算,我们需要为每条路径重新计算dijkstra算法。 在处理大型卷时,批量生成骨架的一个明智方法是对卷进行分块,独立处理分块,并在最后合并生成的骨架片段。然而,这是复杂的"边缘效应"导致的上下文丢失,这使得它不可能期望由相邻块产生的骨架片段的端点对齐。相比之下,很容易连接网格碎片,因为网格碎片边缘的顶点位于给定一个重叠像素的可预测的相同位置。 以前,我们使用50%的重叠来连接相邻的骨架片段,这将大体积骨架化的计算成本增加了8倍。但是,如果我们可以强制骨骼位于边界上的可预测位置,我们可以使用单像素重叠并复制简单的网格连接方法。作为一种(不正确但有用的)直觉,考虑计算每个边界平面上每个连接组件的质心,并将其添加为所需的路径目标。这将保证平面的两边连接在同一个像素上。但是,质心可能不在非凸面外壳的内部,因此我们必须更加复杂,并在形状内部选择一些真实的点。 为此,我们再次重新调整欧氏距离变换的用途,并将其应用于连接组件的六个平面中的每一个,并选择最大值作为强制目标。这对于许多类型的对象都很有效,这些对象与单个平面接触并且具有单个最大值。但是,我们必须处理具有多个最大值的长方体和形状的角。 为了处理与长方体多个侧面接触的形状,我们只需将目标分配给所有连接的组件。如果这在后处理中引入了一个循环,那么我们已经有了在火成岩中处理它的循环移除代码。如果它引入了一些无用的小附件,我们也有代码来处理这个问题。 如果形状具有多个距离变换最大值,则在不需要在可能在不同机器上不同时间运行的空间相邻任务之间通信的情况下选择相同的像素是很重要的。此外,相邻任务上的同一平面会翻转坐标系。一种简单的方法可能是在其中一个坐标帧中选择具有最小x和y(或其他基于坐标的标准)的坐标,但这需要跟踪所有六个平面上的翻转,而且很烦人。取而代之的是,我们使用了一系列基于无坐标拓扑的过滤器,这些过滤器既有趣,又省力,而且选择了一些外观合理的东西。对这种方法的一个有效的批评是,它会在一个完全对称的物体上失败,但这些物体在生物数据中是罕见的。 我们应用一系列过滤器,并根据通过的第一个过滤器选择点: 重要的是,filter 1必须基于标签的形状,以便使凸面外壳的扭结最小化。例如,最初我们只使用2到5个过滤器,但这导致远离大块中心的神经突骨骼突然在区块边界处向区块中心移动。 最初的teasar论文呼吁在步骤6(iii)中使用"滚动无效球"来擦除前景体素。这个球的一个简单实现非常昂贵,因为路径中的每个体素都需要自己的球,而且这些体素中有许多重叠。在某些情况下,可能需要对从根到目标路径上的每个体素无意义地重新评估整个体积。虽然在最坏的情况下可能会出现特殊情况,但在更常见的一般情况下,会花费大量的重复工作。 因此,我们应用了一种使用拓扑线索的算法来执行线性时间内的失效操作。为了实现的简单,我们用立方体代替了球体。函数名 两遍算法如下。给定一个二值图像i、一个骨架s和一组顶点v: 为了使这个模块成为可能,必须对几个经典算法进行特殊调整。 此模块最初设计用于CloudVolume和火成岩。 此软件包中使用的一些teasar修改首先由alex bae演示。 亚历克斯·贝克·戴夫介绍了前体骨架化包和我们在该包中使用的teasar的一些改进。亚历克斯还开发了后处理方法,用于使用50%重叠拼接骨骼。will silversmith将这些技术用于大规模生产,改进了一些基本的算法来同时处理数千个标签,并将它们重写到kimimaro包中。将添加涓流daf,零加权先前探索的路径,并固定边界的算法。forrest collman增加了参数灵活性,并帮助调整daf计算性能。Sven Dorkenwald和Forrest都提供了有益的讨论和反馈。灰尘阈值
最大路径
pdrf_标度
和pdrf_指数
pdrf_scale
和pdrf_exponent
表示惩罚方程的参数,该惩罚方程采用欧几里德距离场(d)并对其进行扩展,以便更接近边界的切割是非常不利的,从而使dijkstra走的路径更居中。pdrf_标度
*(1-d/max(d)pdrf_指数
+(方向梯度<;1.0)。soma_acceptance_threshold
和soma_detection_threshold
soma_acceptance_threshold
是物理半径(例如,以纳米为单位),在此半径之外,我们将图像的连接组件分类为包含soma。距离变换的输出被标签上的空洞所抑制,这些空洞通常是由somata上的分割算法产生的。我们可以填充它们,但是我们使用的孔填充算法很慢,所以我们只希望偶尔应用它。因此,我们设置了一个较低的阈值,即soma_acceptance_threshold
,超过该阈值,我们将填充孔并重新测试soma。soma_invalization_scale
和soma_invalization_const
r(x,y,z)=soma_invalization_scale*dbf(x,y,z)+soma_invalization_const
描述,其中dbf(x,y,z)是距该点形状边界的物理距离。如果做得正确,这可以防止骨骼被绘制到胞体的边界,而是将骨骼主要拉到从细胞体延伸的过程中。修复边框
修复分支
进度
并行
性能提示
对象id
以绕过处理所有其他标签。如果对象id
只包含一个标签,屏蔽操作将运行得更快。cc_safety_factor
<;1来节省高峰内存使用量。pdrf_指数
为2的小功率(例如1、2、4、8、16),以获得较小的加速。fix_branching=false
设置为1.1x到1.5x的中等加速(假设您的teasar参数和数据允许分支)。max_paths
设置为某个合理的级别来保存它们以备以后使用,该级别将在算法检测到至少需要这么多路径后中止并继续下一个标签。动机
为什么是提萨?
teasar衍生算法
一、序言
Ⅱ。骨骼化
fix\u branching=true
pdrf评分
*((1-el/max(el)^pdrf指数
)+dr/max(dr)。scale
*el(v)+const
的无效立方体,并将与这些立方体重叠的il中的任何前景像素转换为背景像素。soma_detection_threshold
..soma_acceptance_threshold
,则转到soma处理模式。soma_invalization嫒u scale
*max(el)+soma_invalization嫒u const
的无效球体,并从其中包含的il中擦除前景体素。这有助于防止在SOMA上绘制出错误的路径。三、定稿
与teasar的偏差
将DAF用于目标,PDRF用于寻路
PDRF = 5000 * (1 - DBF / max(DBF))^16 + DAF
零加权先前路径(
fix_branching=true
)非重叠分块处理(
fix_borders=true
)滚动失效立方体
roll\u invalization\u cube
旨在引起这种尴尬,尽管它似乎并不重要。相关项目
学分
参考文献
推荐PyPI第三方库