Handle image mismatches when only difference is additional singleton axes

[jdtournier]

As discussed in #343

[thijsdhollander]

In this context: I just had a case with an "SH coefficient" file of lmax=0. Depending on whether the 4th singleton dimension is present or not, Math::SH::check will accept or reject it, even though the actual data is of course exactly the same.

[jdtournier]

Yes, that's exactly the kind of daft behaviour I was worried about... But how should this be handled? You'd face the same problem in say Matlab, routinely needing to squeeze() or something to massage the data to fit the expectations of the next stage. We could support some form of what Python calls broadcasting (mrcalc already does), where missing dimensions are treated as having size 1, but I'm not sure that'll be appropriate in all cases.

Feeling a bit short of inspiration on this one. Anyone got any brilliant ideas they'd like to share...?

[Lestropie]

Not sure that we want to accept any 3D volume as input to any command that expects SH volumes, but clearly the case needs to be handled. My instinctive response is that in this specific case the user could run the file back through mrconvert with -axes 0,1,2,3, which would add the singleton dimension; and Math::SH::check() could make the 'suggestion' that if an input 3D image is intended as an lmax=0 volume, this step is required. But that may or may not make sense once I read it back to myself...

[jdtournier]

I've been thinking a bit about this issue, given the recent issues users have been facing, and #552. How about we add a method to the Image class, called say .with_ndim(int), which behaves similarly to .get_image() and .with_direct_io(), in that it returns a new fully formed Image with the number of dimensions fixed up if needed to make the image have the requested dimensionality? The idea would be to use it at Image::open() time, like this:

auto in = Image<type>::open (filename).with_ndim(4);

and optionally with .with_direct_io() tucked in after that. All it would do is return a valid image with the missing dimensions added with size 1 when the dimensionality of the input is less than expected, and throw an Exception if it's higher and the excess dimensions are non-singleton. That should handle case like this OK, right...?

[draffelt]

Sounds like a nice solution to me. I've been caught up with singleton dimensions elsewhere in the past.

[thijsdhollander]

Yep, sounds great. Using this is then the responsibility of the app developer, rather than the end user; which I think makes plain sense.

What about functions like, e.g., Math::SH::check though? Would they need to accept a 3D image? It's nice to be able to run them on a header before pulling out the big guns and .get_image().with_direct_io()-ing the potentially large SH image; so I'd argue they accept a 3D image. In the end, a 3D image is a valid lmax=0 SH image, right?

[jdtournier]

What about functions like, e.g., Math::SH::check though? Would they need to accept a 3D image? It's nice to be able to run them on a header before pulling out the big guns and .get_image().with_direct_io()-ing the potentially large SH image; so I'd argue they accept a 3D image. In the end, a 3D image is a valid lmax=0 SH image, right?

So that might sound nice, but in practice I have a feeling it's not going to be useful. A few reasons:

• most of the SH apps don't currently open a Header before invoking check_SH() - they go for a straight Image::open(). So that would require a bit more modification to the apps than just inserting that .with_ndim() call.
• getting the wrong SH image should be the exception rather than the rule - I reckon users won't mind waiting a little to be told they screwed up on the odd occasion when that happens (if it actually causes a delay at all -see next point).
• most of the time, the image will be accessed via memory mapping - no delay to speak of. Only when loading images like DICOM would there be an explicit load to RAM - unless the app explicitly requests direct IO, which I don't expect will be all that common.
• probably the biggest reason not to allow check_SH() to allow a 3D volume is that it allows the app to proceed with processing, expecting the image to be 4D. So that would either crash the loop or at best result in a dimension mismatch Exception (hopefully the latter...). So it would need to be coupled with a mechanism to ensure the image is also 4D once opened. I'd rather keep things simple...

[thijsdhollander]

most of the time, the image will be accessed via memory mapping - no delay to speak of. Only when loading images like DICOM would there be an explicit load to RAM - unless the app explicitly requests direct IO, which I don't expect will be all that common

Fair enough, I must admit I had not factored this one in specifically.

It's maybe more of a semantic issue: I was more coming from the view that checking for SH contents was just about checking whether the data was technically allowing for an SH interpretation (so basically just checking the number of volumes, and accepting 1, 6, ... as valid cases); and that checking for the 4D-aspect, or even just not be willing to take the lmax=0 case in (because you rely on orientationally varying amplitudes or something), is another separate thing.

But ok, I won't be fussy about it. :-) It's of course also possible just getting the image with_ndim(4) first, running the SH check on that, and then only beyond that point adding in the .with_direct_io().