Compute gradient in physical and voxel coordinate systems

[po09i]

Description

Shimming often requires to compute the gradient of an acquisition. The gradient can be calculated along different axis: the physical x, y, z or the voxel x, y, z. The former will be called Gx, Gy, Gz and the latter Gfreq, Gphase and Gslice for simplicity. Calculating freq, Gphase, Gslice is trivial with np.gradient but calculating Gx, Gy, Gz is not.

This PR introduces a function phys_gradient() that calculates the gradient along the physical coordinates: Gx, Gy, Gz.

Output format

• seperate variables: Gx, Gy, Gz

This PR introduces a function that calculates the gradient along the voxel coordinates: Gaxis1, Gaxis2, Gaxis3 from the physical gradient. phys_to_vox_gradient()

Output format

• seperate variables: Gaxis1, Gaxis2, Gaxis3

Linked issues

Related to #134

[po09i]

I have given a lot of thought about the different gradients and I think some clarification could be useful. From what I have gathered, Gx, Gy, Gz are self explanatory, the change of an input along the physical x, y and z direction. The definition of the gradient along freq, phase and slice axis can be unclear. Here is where i'm not sure of the definition:

Let's take a 1D matrix: data to apply the gradient on = [0, 1, 2] with indexes [0, 1, 2] where the physical coordinates are [0, -2, -4]
Gx:

• ( (1 - 0) / (-2 - 0) ) = -0.5

But what should the gradient be for Gaxis/Gfreq/...

• 1 (does not account for the length of a pixel or physical orientation)
• -1 (does not account for the length of a pixel or but accounts for physical orientation)
• 0.5 (Accounts for the length of a pixel but does not account for physical orientation)
• -0.5

@jcohenadad @evaalonsoortiz

[jcohenadad]

@po09i i think the gradients for the "physical coordinates" (ie scanner) and the "voxel-based coordinates" (ie image) should be the same. At the end, what we care is the effective absolute gradient that the object to image is submitted to. So it should be agnostic to the resolution.

So, i think the example above is missing the actual unit dimensions. To calculate Gaxis, the unit needs to be associated with the matrix, yielding: ( (1 - 0) / (1 - 0) * pix_size ) = 0.5 (with pix_size=2). Now, we have an opposite sign, however, which could probably be taken into account by how the indexation is done (eg: the pix_size should maybe reflect this change in indexation direction?).

[po09i]

@jcohenadad

So, i think the example above is missing the actual unit dimensions.

The physical coordinates are: [0, -2, -4] for [0, 1, 2] indexes. So pix_dim of "-2".

i think the gradients for the "physical coordinates" (ie scanner) and the "voxel-based coordinates" (ie image) should be the same.

I was thinking the opposite actually, since they are in different coordinate systems they would have different gradients in their respective coordinate systems which would result in telling the same effective thing. This is the goal of the functions in this PR; to convert between the two using the affine transformation to link the two coordinate systems. This example shows a negative pix_dim (or rotation of 180 degrees around one axis) but what if there was a 45 degree rotation around the z axis, we would not expect Gx and Gaxis to be the same since Gaxis would be a combination of both Gx and Gy. So for the 180 degrees case, I would expect Gx and Gaxis to be opposite.

I'm leaning towards an answer for the example above to be either 1 or 0.5 actually. Either as if the indexes were the coordinates of the image or that plus accounting for pix_dim but using the axis orientation. I'm not confident enough in my reasoning though it might depends on what it means when Siemens expects Gaxis.

[jcohenadad]

I was thinking the opposite actually, since they are in different coordinate systems they would have different gradients in their respective coordinate systems which would result in telling the same effective thing. This is the goal of the functions in this PR; to convert between the two using the affine transformation to link the two coordinate systems. This example shows a negative pix_dim (or rotation of 180 degrees around one axis) but what if there was a 45 degree rotation around the z axis, we would not expect Gx and Gaxis to be the same since Gaxis would be a combination of both Gx and Gy. So for the 180 degrees case, I would expect Gx and Gaxis to be opposite.

I agree with you about the orientation (ie Gslice and Gz might not be the same), however the norm of the applied gradient should be the same. If they were not the same, it would mean that depending on the slice thickness of the target image, the applied Gslice would vary (whereas the effective static B0 gradient is the same).

[jcohenadad]

I'm approving without extensively digging into the code as it would take me some time to double-check everything. I'm glad you added some tests, which offer a first layer of confidence. I think we should also add another layer of testing with real data acquired in various orientations. Great job 👍