Questions regarding the computation of the path lenght

[ReubenDo]

Hi Guotai,

I have few questions regarding the way you use the spacial information in your code:

GeodisTK/cpp/geodesic_distance_3d.cpp

Line 287 in 7907b38

float delta_d = local_dis_b[i] / speed;

I am not sure that this factor is correct, I think that the factor should be proportional to 1/local_dis_b[i].

Let's be given an image I with a spacing [a, b, c] and a point S used as seed. S is located at [x,y,z] in the voxel space and at S=O + [xa, yb, z*c] where O is the origin in the scanner space.

Let's consider that we want to know the length of the straight path between S and the voxel just above in the z axis, i.e. T = O + [xa, yb, (z+1)*c].

Then, I would say that the length of the path is:

where the unit vector u = Γ ′(s)/||Γ ′(s)|| is tangent to the direction of the path.

If you agree with that, that would also change the way the normalisation is performed for the other path direction.
For example if the curve has a direction (1,1,1), this would lead to a factor

instead of:

Please let me know what you think.
Thank you!

Regards,
Reuben

[taigw]

Hi Reuben,
In order to make Euclidean distance and Geodesic distance can be formulated in a single equation, this code uses a general idea that the total length of a path is the sum of all the segments of that path in the physical space. The length of each segment is defined by its spatial length divided by propagating speed, which is similar to the fast marching idea. In a special case, if we define the speed as a constant (speed = 1), the output is exactly the spatial length of that path. In the geodesic case, the speed is 1/gradient.
For your question, I think the result in your first example should be multiplied by c, so that the value is measured in the physical space, which is independent of the resolution of the image. For example, if we consider three points in the x axis, A(x=0 mm), B(x=1 mm), C(x=2mm), and the intensity values of them are I_A = 0, I_B = 100, I_C = 200. We want to calculate the geodesic distance between A and C. If we have two images for the same physical objects: image1 has a resolution of r1=1mm, and image2 has a resolution of r2=2mm. Then d(A,C) in the first image should be d(A,B) + d(B,C) = 100/r1 * 1 + 100/r1 * 1 = 200, and d(A,C) in the second image should be 200/r2*2 = 200. So the result in the physical space is intendent of the imaging process.
Guotai

[ReubenDo]

Hi Guotai,

Thanks for your answer. Great example.

I don't understand why d(A,C,Image1) should be equal to d(A,C,Image2) when the image gradient is used. For example, we could imagine that I_B = 10000 and I don't see why the cumulative sum of the gradient should be equal in both cases.

In contrast, if we consider that A and C are directly connected in both case, then the length of the two direct paths should be equal.
For image 1, the image gradient at C is (I_C-I_A)/(2.r1) (factor 2 because gap of 2 voxels)
For image 2, the image gradient at C is (I_C-I_A)/r2.
And we have (I_C-I_A)/(2.r1) = (I_C-I_A)/r2 with r1=1 and r2=2.

To test that, I use SimpleITK to compute image gradients using the spacing information. Then I computed the gradient using numpy that uses the first-order finite difference:

The example below confirms that the factor used in SimpleITK to compute the gradient is 1/spacing[i]:

import numpy as np
import SimpleITK as sitk 

input_name = "data/img3d.nii.gz"
img = sitk.ReadImage(input_name)
I = sitk.GetArrayFromImage(img).transpose()

# Compute Image gradient using SITK
image_grad= sitk.GradientImageFilter().Execute(img)
grad_sitk = sitk.GetArrayFromImage(image_grad).transpose()

# Manuel computation of the gradient using numpy:
spacing_raw = img.GetSpacing()

grad_np = []
for i in range(3):
    factor = 1/spacing_raw[i]
    grad_np.append(factor*np.gradient(I, axis=i, edge_order=1))
grad_np = np.stack(grad_np,0)


# Remove borders
grad_np_noborder = grad_np[:,2:-2,2:-2,2:-2]
grad_sitk_noborder = grad_sitk[:,2:-2,2:-2,2:-2]
diff = abs(grad_np_noborder-grad_sitk_noborder)
print(f"Gap between techniques: f{diff.mean()}")

Gap between techniques: f8.152106643539441e-07

Reuben

[taigw]

Thanks Reuben,
I think you are right. The previous code ignored the spacing information when calculating the gradient. This can be fixed by updating line 286 as follows:
float speed = (1.0 - lambda) + lambda* local_dis_b[i]/(l2dis + 1e-5);
where local_dis_b[i]/(l2dis + 1e-5) is the inverse of gradient. The forward pass can be updated in the same way.
Guotai