add set functions to GeometricalInfo

[KrisThielemans]

It would be intuitive to be able to manipulate a GeometricalInfo, and then create a new object from it (#824 ). You might want to change the voxel-size, or extent/bounding box, or even offset.

[KrisThielemans]

before we go ahead and implement it, let's discuss an interface for it here.

[DANAJK]

In some of my other work, I use a structure to keep track of geometry. This has fields that parallel the DICOM geometry description (3D patient coordinate of centre of top left voxel, pixel sizes (mm) and the 3D direction of the rowa and columns). I then have a function that will update these parameters for in-plane 'transformations' such as cropping, padding, voxel size change, rotation and flipping. Having this is quite liberating because you are free to manipulate images whilst knowing that the geometry information remains correct.

[KrisThielemans]

@DANAJK, exactly. That's what 'GeometricalInfo' is (except that it is a class and not a structure), C++, Python, Matlab

However, at present we cannot manipulate it. It's initialised for an ImageData when reading, and we can get info from it.

I guess you're saying you want things like new_geo_info=GeometricalInfo.crop(..)). It'd be great to agree on a good interface for this (names, parameters).

[DANAJK]

For many applications, you will probably want to update an existing Geometrical Info, rather than make a new copy.

The names I used for the methods were: 'rot', 'flip', 'circshift', 'remove_oversampling', 'zeropad' and 'imresize'. I update the image data and the geometrical data within the same function - this keeps the code in one physical file.

[ashgillman]

I'm not quite sure when you might want to change voxel size (you should be using the resampling facility for this, right?).

Changing the extent and bounding box can't be done without also modifying the image container, and could also be done with the resampling operator.

I think the main reason you'd want to do this is to apply a transform to an image without resampling (an operation we called "reorient" in the Mirorr project).

I kind of like having the GeometricalInfo class as static - I'd rather see an: image_data.reorient(transform) # replaces image_data's _geom_info_sptr internally function. I see there is already such a function, but I'm struggling to work out what it does - it might need a bit more explicit documentation.

[DANAJK]

You need to change voxel size when comparing two reconstructions that used different underlying voxel sizes, e.g. in MR a high resolution T2 image and a diffusion image. For PET/MR, this could be to compare the PET and MR reconstructions.

I agree this needs a resampling. I don't mind where the code that updates the geometrical info goes, so long as it is documented and tightly coupled to the code that changes the data (so that you can't easily update one and not the other.)

For resizing, updating the geom info is trivial if you are resampling image A to the geometry of image B because the geom info can then be copied from B. It is more complicated if you just call a generic resampling function (e.g. imresize in MATLAB) because the 3D patient position of the top left voxel will differ in the resampled data - it is not enough to just update the voxel size.

[KrisThielemans]

I personally prefer to make resampling explicit. It also allows you to specify interpolation, boundary conditions etc. Of course, our current interface is too verbose for many, but here it is in full generality as an example

resampler = sirf.Reg.NiftyResample()
resampler.set_reference_image(my_template_target)
resampler.set_floating_image(PET_mrac)
resampler.set_interpolation_type_to_linear()
def=sirf.Reg.NiftiImageData3DDeformation('transf%d.nii' % g)
resampler.add_transformation(def)
resampled = resampler.forward (PET_mrac)

If you would use an identity transformation, this would keep geometry in the same place. I guess at present we have to create an explicit identity transformation, but that'd be easy to solve (I guess it could even be made the default). And I'm always confused about the reference/floating image terminology, but I guess you can see that this could be easily shortened to

resampler = sirf.Reg.NiftyResample(target_geometrical_info=my_template.get_geometrical_info())
resampled = resampler.forward (PET_mrac)

and then of course we can get rid of the explicit call to get_geometrical_info.

So, I was thinking that if we can manipulate the geom_info, then we can rotate etc, as in

my_target_geo = PET_mrac.get_geometrical_info().rotate(whatever_args)
resampler = sirf.Reg.NiftyResample(target_geometrical_info=my_target_geo)
resampled = resampler.forward (PET_mrac)

and yes, this can be shortened as well, but I hope you get the idea.

[DANAJK]

OK. So after this, would resampled.get_geometrical_info() correctly return the updated geometrical info?

[KrisThielemans]

yes it does that already now.

[KrisThielemans]

If desired, we could add the rotate et al methods to ImageData I guess, with an optional specification of the resampler. It would then do the above lines. It would lead to more code duplication at the class level (but possibly less for the user).

[KrisThielemans]

Note that there are complications here about the type of ImageData. What I've lost with the target_geometrical_info option is the type of the desired output image. I guess can be folded in as well.

Another caveat is that with named arguments in Python, things are rather clear, but neither C++ nor Matlab have that, so there things are messier.

[DANAJK]

The example below would be fairly intuitive, but it needs the rotate method to have an implicit understanding of how NiftyResample is going to interpret origins, rotation signs etc?

my_target_geo = PET_mrac.get_geometrical_info().rotate(whatever_args)
resampler = sirf.Reg.NiftyResample(target_geometrical_info=my_target_geo)
resampled = resampler.forward (PET_mrac)

[ashgillman]

I don't think it's necessary to pass the geom_info to the resampled in your shortened example. Why not:

resampler = sirf.Reg.NiftyResample(reference_image=my_template)
resampled = resampler.forward (PET_mrac)

I think this would be a lot easier to implement (I think Nifty internally would need an image) and solves the data type problem. And seems intuitive to anyone who has used registration tools before.

I think the image_data.reorient(transform) proposal above would also accommodate what you are suggesting.

We could add convenience operations like:

rotation_transform = Transform(rot_x_degrees=45, trans_z_mm=-50)

Would require some thought as to the Transform interface, but it allows less working parts (just a transformer for the image)

[ashgillman]

So to rotate an image 45 degrees about X without resampling:

image.reorient(Transform(rot_x_degrees=30))

[KrisThielemans]

it needs the rotate method to have an implicit understanding of how NiftyResample is going to interpret origins, rotation signs etc?

Everything in SIRF related to geom_info is supposed to work in LPS. We can represent rotations in different way as we know of course. It'd probably explicitly specify rotation centre, as otherwise people will get confused.

I don't think it's necessary to pass the geom_info to the resampled in your shortened example. Why not:

resampler = sirf.Reg.NiftyResample(reference_image=my_template)
resampled = resampler.forward (PET_mrac)

yes, of course. This is essentially our current way of doing this. However, this has the eternal problem that you need to have a template-image. With this issue, I was aiming to get round that problem by finding good ways to create that. Once we've resolved that, and after #824, we could say

my_template=NiftiImageData3D(my_target_geo) # or STIRImageData or ....
resampler = sirf.Reg.NiftyResample(reference_image=my_template)
...

or we can decide to hide this in the NiftyResample constructor options.

Of course, there's probably a million ways to do this, but I think we should not provide too many paths, as that leads to confusion as well. (IMHO it's better to have a consistent (and well-documented) interface, than one with that allows one-liners, but that are different for every case).

We could add convenience operations like:

rotation_transform = Transform(rot_x_degrees=45, trans_z_mm=-50)

let's not introduce x,y,z terminology please but stick to LPS.

However, that made me think that we have AffineTransformation already. This has roughly what you have above.

So, if we can apply an AffineTransform on a GeometricalInfo, and solve #824, are we nearly there?

[ashgillman]

Perfect, sorry I'd have looked up the existing implementation naming but I'm just on mobile.

I think so.. we didn't quite discuss cropping, which works with resampler but I guess shouldn't really need building with Riftyreg..