Apply global rotation to (axis) geometries

[kohr-h]

When reconstructing arbitrary slices, one can trick ODL into doing that by reconstructing in a space with shape (n0, n1, 1) and rotating the whole geometry with the matrix rotating the normal vector of the slice to the original geometry axis. Of course, the same can be done in 2D, too. To do it right, one must rotate the following geometry vectors:

• axis (if present)
• src_to_detector_init or det_init_pos
• det_init_axes or det_init_axis

I think a helper function or method would be good to get the rotations right -- what I'm not so sure about is whether this is functionality we need in the core. Opinions?

[adler-j]

Frankly I'm not sure if this should be in core or not. I feel that there is a quite limited applicability while complicating the interface with questions like how does this interact with other parameters.

Why not simply change the angle in the projections? Are you using a tilted geometry or something?

[kohr-h]

I feel that there is a quite limited applicability

True. But there's more of that to come.

complicating the interface with questions like how does this interact with other parameters.

Interface doesn't change, everything goes via init parameters. You just make a new geometry from the old one, rotating the input vectors in the correct way.

[adler-j]

Well my second question still holds. How is this different from simply giving an angle_partition that is slightly offset? I don't see the use case right now.

[ozanoktem]

@kohr-h , is the need for having a utility to generate a parallel beam geometry along an arbitrary rotation axis related to the tilt-axis orientation in single-axis tilting in ET?

[kohr-h]

@adler-j It's totally different. If you offset the angle partition, you get the same geometry starting and ending at different angles. If you have an axis-oriented geometry, that's an intrinsic rotation around the axis. What I'm talking about is rotating the axis, the source curve and the detector curve all simultaneously, thus applying an extrinsic rotation. Basically it's a transformation of the "world" system, but we keep sitting in it so for us the geometry rotates. @ozanoktem Not quite, in that case it's only the tilt axis that's rotated while the detector keeps its orientation.

[adler-j]

I still don't see that you are aiming at here, could you make a figure or something?

And regarding the design, I guess your idea is to make a helper function for this? What would you call it and what would the parameters be? perhaps that would help.

[kohr-h]

Sure, check out the code in this example. The workflow is as follows:

• Define global geometry and ray trafo
• Generate some projection data
• Filter the data
• Define a slice by giving a normal vector and a shift
• Create a reconstruction space with shape (n0, n1, 1) and shifted corners. This can be understood as the "world" rotation from the original slice normal to the new normal (0, 0, 1).
• Create a geometry that is rotated in the same way, thus keeping the relative orientation between slice and geometry constant.
• Create a ray trafo using the "slice space" and the rotated geometry and apply its adjoint to the filtered data.

Result: FBP reconstruction in the rotated and shifted slice only.

Some images:

Standard orientation (x-y plane)

x-z plane

y-z plane

Now for something skew - normal vector (1, 1, 0)

Shifting along the new axis lets us go "up and down"

Shifts not along the new axis are lateral

[kohr-h]

Note that the last 3 points contain quite some boilerplate code which I feel should be put into some helper function. There are rotations all over the place, and it's easy to get it wrong (matrix or transpose? what to rotate? rotate then shift or the other way around?) so I would like to make it easier.

This is related to #783 and #152. Most of the time here is spent swapping the array axes and copying the data to GPU. The actual BP only takes 70 ms on my machine, that would pretty much give you 15 FPS "live reco".

[adler-j]

How is what you did different from this code? I feel like I am missing something still. The only difference is that you give the shift in the new coordinate system, while I give it in the old system, and you seem to be applying a different rotation than I do, but that is to be expected since determining a rotation given a "from" vector and a "to" vector is under-determined anyway.

Example of the code running your second to last example:

[kohr-h]

How is what you did different from this code? I feel like I am missing something still. The only difference is that you give the shift in the new coordinate system, while I give it in the old system, and you seem to be applying a different rotation than I do,

That's the point. The initializer of the geometry chooses the detector coordinate system arbitrarily and not consistent with the rotation of the axis, i.e. the detector orientation is not the rotated version of the old one, except for special cases. The same applies to src_to_det_init, it's just some arbitrary perpendicular_vector to the chosen axis,

but that is to be expected since determining a rotation given a "from" vector and a "to" vector is under-determined anyway.

which is also under-determined.

So in summary, I apply the rotation consistently, you rotate the axis and rely on defaults otherwise.

[kohr-h]

And I'm thinking of adding an intrinsic rotation that aligns the slice axes with the "image axes". Then it's a unique rotation.

[adler-j]

I may be missing some definitions here, but in general picking a rotation given a "from" and a "to" is under determined up to a (one dimensional) rotation parameter, so the only difference between our codes would be how that rotation is picked.

You seem to have another choice than what ODL does, but in general ANY choice will be a "bad choice" due to the hairy ball theorem, i.e. you always get a singularity somewhere, where a minimal change in the "to" axis drastically changes the rotation. For you, that happens when "from" and "to" are close to parallel/anti-parallel, while in ODL this happens when either axis is close to +- [0, 0, 1].

Anyway, what I'm getting at is that if the difference between your approach and the one I outlined is down to the specific choice of that free rotation parameter, wouldn't it be better to focus on that instead of adding a whole machinery for this? I guess we could simply add a text parameter that specifies how the other axes should be selected given an axis. Then you could add your idea in addition to the current one (or perhaps simply change the current implementation)

[kohr-h]

I may be missing some definitions here, but in general picking a rotation given a "from" and a "to" is under determined up to a (one dimensional) rotation parameter, so the only difference between our codes would be how that rotation is picked.

I still think that's not the case, rather that your version rotates the axis and then picks the default "derived vectors" while I take the derived vectors in the old system and rotate them all simultaneously.

You seem to have another choice than what ODL does, but in general ANY choice will be a "bad choice" due to the hairy ball theorem, i.e. you always get a singularity somewhere, where a minimal change in the "to" axis drastically changes the rotation. For you, that happens when "from" and "to" are close to parallel/anti-parallel, while in ODL this happens when either axis is close to +- [0, 0, 1].

See above, that's not the problem I try to solve. I want to rotate everything the same way, no more, no less.

Anyway, what I'm getting at is that if the difference between your approach and the one I outlined is down to the specific choice of that free rotation parameter, wouldn't it be better to focus on that instead of adding a whole machinery for this? I guess we could simply add a text parameter that specifies how the other axes should be selected given an axis. Then you could add your idea in addition to the current one (or perhaps simply change the current implementation)

So what would be an option is to internally define defaults for all vectors based on the axis=[0, 0, 1] case, that is src_to_det_init=[1, 0, 0] and det_init_axes=([0, 1, 0], [0, 0, 1]) and then rotating everything with the same matrix for other choices of axis, which is what I'm doing externally right now. For me that would be a good solution since my code would reduce to making a new geometry just with a rotated axis.

[adler-j]

I still think that's not the case, rather that your version rotates the axis and then picks the default "derived vectors" while I take the derived vectors in the old system and rotate them all simultaneously.

I fully agree that there is a difference in how you think about it, but looking at your code it all boils down to the "rotation_matrix" function, and in there you make a particular choice. The "ODL" choice would satisfy the same conditions i.e. transfer "from_vec" to "to_vec", but would give a different matrix.

So sure, the casuality arrows point in different directions, but the correlation is still the same.

See above, that's not the problem I try to solve. I want to rotate everything the same way, no more, no less.

My code also rotates everything the same way, so there has to be some other constraint that you add that I am missing, and I think this is where the confusion is. Perhaps it is related to your comment "And I'm thinking of adding an intrinsic rotation that aligns the slice axes with the "image axes". Then it's a unique rotation.", but I don't think the current code does that.

So what would be an option is to internally define defaults for all vectors based on the axis=[0, 0, 1] case, that is src_to_det_init=[1, 0, 0] and det_init_axes=([0, 1, 0], [0, 0, 1]) and then rotating everything with the same matrix for other choices of axis, which is what I'm doing externally right now. For me that would be a good solution since my code would reduce to making a new geometry just with a rotated axis.

Morally, that is what the current code does. Now, the implementation is not the best and actually writing this using a rotation matrix would likely be much preferable to the current way of writing it, I agree with that. And in doing that, we would get this one free parameter that I've been talking about here, and we could give users an option in how to select that.

[kohr-h]

I think what it all boils down to is that the perpendicular_vector function applies a rotation (if it can always written as one) that is not a rotation from [0, 0, 1] to axis. That's the constraint I add, and we could choose to apply that in general.

[adler-j]

I think what it all boils down to is that the perpendicular_vector function applies a rotation (if it can always written as one) that is not a rotation from [0, 0, 1] to axis. That's the constraint I add, and we could choose to apply that in general.

I actually think that it does, since the free variable is "orthogonal" to "being a rotation from [0, 0, 1] to axis" in some sense.

Anyway an obvious way to do this would be to stop using the det_init_axes parameters as much internally and instead move the logic into finding a rotation matrix and pre-multiplying with it here and in other related places.

We could even let users provide the full rotation matrix themselves as an option to providing axis.

Finally, the current "det_init_axes" should then ONLY be used for when users want detectors with non-rectangular pixels, or when the detector should not point straight towards the source.

[kohr-h]

I actually think that it does, since the free variable is "orthogonal" to "being a rotation from [0, 0, 1] to axis" in some sense.

>>> unit_vecs = np.eye(3)
>>> perp_vec_matrix = np.zeros((3, 3))
>>> for i in range(3):
...     for j in range(3):
...         perp_vec_matrix[i, j] = perpendicular_vector(unit_vecs[j]).dot(
...             unit_vecs[i])
...
>>> print(perp_vec_matrix)
[[ 0. -1.  1.]
 [ 1.  0.  0.]
 [ 0.  0.  0.]]

Doesn't look like a rotation matrix to me.

Anyway an obvious way to do this would be to stop using the det_init_axes parameters as much internally and instead move the logic into finding a rotation matrix and pre-multiplying with it here and in other related places.

We could even let users provide the full rotation matrix themselves as an option to providing axis.

Finally, the current "det_init_axes" should then ONLY be used for when users want detectors with non-rectangular pixels, or when the detector should not point straight towards the source.

Yes, I think that's a decent suggestion. Perhaps good to let it settle a bit over Christmas :-)

[adler-j]

Doesn't look like a rotation matrix to me.

Sure several calls to this will not give a rotation matrix, in fact it is mathematically impossible to write a function that takes a vector and returns another vector that is perpendicular to the initial vector in a continuous way, and if I'm not wrong that rules out a matrix since a matrix would be linear and thus continuous.

I could get the same problem if I called your code with from_vec == to_vec, for example.

However, the function is only called once for each geometry, and thus that does not become a problem.

Yes, I think that's a decent suggestion. Perhaps good to let it settle a bit over Christmas :-)

Yupp. I think this is best settled with an overall look at geometries, they are not top notch right now.

[kohr-h]

That's true, and it's not really the point of the discussion. However, the resulting geometrical layout of the acquisition system is different for the different approaches. Here is some example code that simulates both approaches. The output is:

------ Default ------
axis: [0, 0, 1]
src_to_det_init: [ 1.  0.  0.]
det_init_axes: [array([ 0.,  1.,  0.]), array([ 1.,  0.,  0.])]
------ New axis ------
axis: [0, 1, 0]
src_to_det_init: [-1.  0.  0.]
det_init_axes: [array([ 0., -0.,  1.]), array([-1.,  0.,  0.])]
------ Global rotation ------
axis: [0, 1, 0]
inverse rotation of axis: [  0.00000000e+00   6.12323400e-17   1.00000000e+00]
src_to_det_init: [ 1.  0.  0.]
det_init_axes: [array([  0.00000000e+00,   6.12323400e-17,   1.00000000e+00]), array([ 1.,  0.,  0.])]

And it's not really surprising that the results differ. How should the perpendicular_vector function know what rotation it should apply? It makes some guess that is guaranteed to work, but it can be consistent with a given external rotation only by chance.

So that issue could simply be solved by using a rotation matrix internally instead of the coordinate system(s) set in the initializer.

[adler-j]

However, the resulting geometrical layout of the acquisition system is different for the different approaches.

I agree with that, but my point is that sadly any approach we pick will be arbitrary so if we do something about this it would feel better if they were treated equally (i.e. the current is not magically more important than the one you suggested)

And it's not really surprising that the results differ. How should the perpendicular_vector function know what rotation it should apply? It makes some guess that is guaranteed to work, but it can be consistent with a given external rotation only by chance.

Well they differ by design, but they are both still "guaranteed" to give "rotation matrices" (noticed now that the current implementation does not preserve handedness though), in your case

------ New axis ------
[[-1.  0.  0.]
 [ 0.  0.  1.]
 [ 0.  1.  0.]]
------ Global rotation ------
[[ 1.  0.  0.]
 [ 0.  0.  1.]
 [ 0.  1.  0.]]

So that issue could simply be solved by using a rotation matrix internally instead of the coordinate system(s) set in the initializer.

Agree wholeheartedly.

[adler-j]

Fixed by #968