只有当位置发生变化时，raster_geometry.sphere才会提供不同的卷

Question

我有一组中点(从一组扫描的标记肿瘤到512x512x512体素)。我只是想在每个中点创建一个半径为10个体素的球体。为此，我使用raster_geometry.sphere并将位置定义为中点/ 512 (即位置=中点/512，raster_geometry.sphere(512，10，positon) )

但是，当我这样做的时候，球体的体积会发生变化。关于raster_geometry，有什么我不明白的吗？

from logzero import logger as logging
from lama import common
import numpy as np
import SimpleITK as sitk
import raster_geometry as rg

for i, img_path in enumerate(scan_paths):
    logging.info(img_path)
    logging.info(tumour_paths[i])

    m_loader = common.LoadImage(tumour_paths[i])
    mask = m_loader.img

    m_array = sitk.GetArrayFromImage(mask)

    s = ndimage.find_objects(m_array)[-1]

    midpoint = [(np.mean([s[0].start, s[0].stop]))/512,
                (np.mean([s[1].start, s[1].stop]))/512,
                (np.mean([s[2].start, s[2].stop]))/512]

    print("Original Midpoint", [i*512 for i in midpoint])

    print("Modified midpoint", midpoint)

    arr = rg.sphere(512, 10, midpoint).astype(np.int_)

    print(np.count_nonzero(arr))
    print(np.sum(arr))

示例输出如下：

[I 220525 17:37:58 radiomics_normaliser:86] E:\220204_BQ_dataset\220521_BQ_norm\imgs\200721_MPTLVo3_GFSeeds_4T1R_4T1R_D7_C1_002.nrrd
[I 220525 17:37:58 radiomics_normaliser:87] E:\220204_BQ_dataset\220521_BQ_norm\tumour_respaced\200721_MPTLVo3_GFSeeds_4T1R_4T1R_D7_C1_002.nrrd
Original Midpoint [260.5, 252.5, 162.0]
Modified midpoint [0.5087890625, 0.4931640625, 0.31640625]
4160
4160
[I 220525 17:38:14 radiomics_normaliser:86] E:\220204_BQ_dataset\220521_BQ_norm\imgs\200721_MPTLVo3_GFseeds_4T1R_4T1R_D7_C1_003.nrrd
[I 220525 17:38:14 radiomics_normaliser:87] E:\220204_BQ_dataset\220521_BQ_norm\tumour_respaced\200721_MPTLVo3_GFseeds_4T1R_4T1R_D7_C1_003.nrrd
Original Midpoint [219.5, 234.5, 165.5]
Modified midpoint [0.4287109375, 0.4580078125, 0.3232421875]
4165
4165

Answer 1

这是光栅化的一个基本问题。

让我们试着用ASCII图片来说明这一点。

假设我们想表示一条长度为2的线，从0开始，在光栅中使用4单元：

        0       1       2       3       4
raster: |_._._._|_._._._|_._._._|_._._._|
  line:  - - - - - - - - 

如果单元格只能具有布尔值，则当我们计算"True“单元格的数量并将其用作行长度的代理时，就会得到2。好的。

如果行从0.5而不是0开始，会发生什么？

从图形上看，这是：

        0       1       2       3       4
raster: |_._._._|_._._._|_._._._|_._._._|
  line:      - - - - - - - -

现在这条线跨越3个细胞。如果我们将半填充的单元格计数为"True"，则计数方法将给出一行长度为3的值。相反，如果我们将半填充的细胞计算为"False"，则得到1。

那我们怎么才能解决这个问题？在这种情况下，我们可以说，我们把一个单元格算作“真”，一个算成“假”，然后我们就收支相抵。

然而，这并不理想。如果我们的线的长度是2.25开始呢？无论我们在计算“真”细胞方面有多聪明，我们都不会得到正确的结果。

一种更有效的方法是不仅存储单元格是否填充(部分)，而且存储填充量(从0到1)，即使用小数加权而不是布尔值。

使用上述示例中的2.25行长度和其余部分，您将得到：

        0       1       2       3       4
raster: |_._._._|_._._._|_._._._|_._._._|
  line:  - - - - - - - - -

所以，我们存储在细胞中：[1, 1, 0.25, 0]。

如果行现在从0.5开始，这也不是问题，因为单元格内的值会告诉我们在每个单元格中有多少行：

        0       1       2       3       4
raster: |_._._._|_._._._|_._._._|_._._._|
  line:      - - - - - - - - -

准确地说是：[0.5, 1, 0.75, 0]，它又一次概括为2.25。

对于行来说，产生这样的估计是相对容易的。

对于复杂的形状，这就不那么简单了，并且在raster_geometry中提供了一种尝试性的方法，并且只有当smoothing参数设置为True时才激活它。

另一种方法是用更小的单元取代较大的单元，即更高的分辨率。这样，就可以以内存为代价来实现精度。

在ASCII图片中，如果每个单元格被分成四个，那么2.25大小就变成2.25 * 4 = 9

        0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6
raster: |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|
  line:  - - - - - - - - -

对于线性问题，高分辨率的内存“代价”是：size(bool) / size(float)。

对于64位机器上的NumPy，这是8 / 64 = 8，所以更高的分辨率将导致更粗的长度估计，高达单元大小的1/8。

对于二维，我们需要取平方根，对于三维，对于立方根，对于n，对于n-根，我们需要取平方根。例如，在3D中，相同的内存将产生1/2的单元格大小精度。

当然，如果分辨率更高，就会得到额外的信息，说明额外的重量在哪里，所以上面的讨论只适用于测量物体的尺寸(以及相关的量，如面积、体积等)。