MRtrix3 / mrtrix3

MRtrix3 provides a set of tools to perform various advanced diffusion MRI analyses, including constrained spherical deconvolution (CSD), probabilistic tractography, track-density imaging, and apparent fibre density
http://www.mrtrix.org
Mozilla Public License 2.0
291 stars 179 forks source link

Image::Buffer*: Interface & naming #211

Closed Lestropie closed 9 years ago

Lestropie commented 9 years ago

Thijs is yelling in my ear and won't let me go home until I write this up...

So he ran into a few understanding hiccups on figuring out MRtrix3's image access classes, and got me wondering about whether things could be improved a little bit. Let's enumerate these for the sake of responses:

  1. Naming of Image::Buffer; I'm not a big fan of having Image::Buffer then Image::Buffer*, as it appears as though the former is a base class of the latter, which isn't the case; it's a different mechanism for providing access to image data. Maybe Image::Buffer should in fact be Image::BufferDisk; that way it's explicit both that it provides access to image files, and that it lives alongside e.g. BufferScratch rather than above/below?
  2. Naming of Image::Buffer (again): a 'buffer' could be thought of as a block of data; in our case, image intensities. But our Image::Buffer* also contains header information. So maybe these should instead be called Image::Image*?
  3. I'm sure I've discussed this before: Having a single class that is not templated, but will operate on any type of underlying Image::Buffer* and any data type. I still don't see any way to make it possible (at least not without an absurd amount of code branching), but he's pleading.
jdtournier commented 9 years ago

Man, you guys are going to kill me... This chart springs to mind...

OK, here goes:

Maybe Image::Buffer should in fact be Image::BufferDisk

Possibly, although Image::BufferPreload is also used to access images on disk... A better name would be Image::BufferRaw, which at least hints at the fact that this buffer gives you access to the data as-is, in contrast to Image::BufferPreload.

So maybe these should instead be called Image::Image*?

:) This is a joke, right? I can't see how that's any less confusing...

In a previous version, this used to be called Image::Object, which was probably not a great deal more informative - but at least it avoided all the confusion that would ensue for developers and compiler alike if classes have the same name as their enclosing namespace... Not against a different name, but you'd have to come up with a more convincing one.

The other option is to keep the buffer and header separate, with the buffer only keeping a reference to its parent header. Personally, I think this would be a recipe for disaster: you really don't want to allow the possibility of the header being modified after the buffer is instantiated - this is why the Buffer classes inherit from ConstHeader rather than Header. Besides, I don't think it would change actual use of the API in any way, this would be purely a conceptual difference, with very little practical consequences for the way the code would be written. If you were to keep the Header and Buffers entirely separate, you'd need to provide the relevant information to the Voxel class too (so it would need to be constructed with both the Header and the Buffer), and would hence behave differently with the BufferScratch class, which does not require a full header, it only needs an InfoType interface. This would make a lot of the generic algorithms that we have very awkward... Otherwise, you could make all the Buffer class only provide a limited InfoType interface, constructed from the Header where relevant, so that there is enough information in there for the Voxel class to operate, but by that stage you might as well include the full header...

We actually thought about this name quite a bit, and I really do think Buffer is the best term here. This class provides a view into the image data as a linear array - that's its primary purpose. To be able to access this information meaningfully though, you need to have all the information related to the image dimensions and strides, which is what InfoType provides. And since InfoType is essentially a minimal header anyway, I don't really see a good argument for not including the full header while you're at it...

Having a single class that is not templated, but will operate on any type of underlying Image::Buffer* and any data type. I still don't see any way to make it possible (at least not without an absurd amount of code branching)

Precisely. You can't do this without virtual methods or code branching, and if you do that you can't avoid a significant loss of performance. For instance, when I modified tckgen to use an Image::BufferPreload rather than an Image::Buffer, this improved performance by a factor of ~3 (at least on my laptop, when I was developing this class). This is because you don't get the cost of a full-blown function call per image intensity read/write (and there's a lot of them in there with all the interpolation that goes on, particularly for iFOD2), and the read/writes themselves don't involve any type conversion or scaling.

What I would much rather do is figure out a cleaner API to avoid the need to write template type specifications everywhere, which can be done using a variety of techniques - especially in C++11 with the auto and decltype keywords. What also helps is writing helper functions that return the fully-formed type, so you can declare it with auto - much like the Buffer classes' new .voxel() call.

Anyway, open to ideas here, but only within reason...

draffelt commented 9 years ago

OK, here goes. I'm sure Donald will have good points against all these too. I know we discussed it at length, however that was a while ago and so I've probably forgotten a lot of the issues. Firstly, I don't feel too strongly about any of this, I'm pretty happy with how things are now.

I agree that Image::Buffer may be mistaken as a base class for Image:BufferScratch etc. I do like the idea of making it explicit that it is mapped to disk (something that most people are not used to dealing with). If we changed it to something like Image::BufferDisk, is there any reason why we can't include preload as an option in the constructor (either a bool or by passing the strides), and get rid of the BufferPreload class?

With respect to #2, I can see what Rob is getting at. Strictly speaking a buffer is really an image, since it also contains the header info, and it can now be constructed directly from a file name (without creating a header first). The problem is the namespace issue...I agree that Image::Image is a bit weird. Although Image::ImageDisk, Image::ImageScratch, Image::ImageSparse is less weird.

Alternatively we could remove the image namespace (like ITK's image class is just itk::Image). So MR::Image::Buffer would be MR::Image (or MR::ImageDisk, MR::ImageScratch, and MR::ImageSparse), and then all the filters would be MR::Filter::Dilate etc. Since we (currently) don't have filters, interpolation, or registration that work on anything but images (for example meshes), then it's probably not such a big deal removing the Image namespace. Although things like MR::Handler::Base might need the namespace changed to MR::ImageHandler::Base.

When I first started using MRtrix I also found it weird/redundant that I could access the same header info from either a header, buffer or voxel object (although now I just find it useful). At the time, we had to create a header first before creating a buffer, then we would pass voxels around to various functions. Now that we can directly create and pass around buffers (images!), then create voxels whenever needed using auto, does a voxel really need to be able to return buffer/image attributes such as ndim, dim etc?

jdtournier commented 9 years ago

Ah crap. This has now turned into a full blown mutiny...

OK, this is going to be a long one - better grab yourselves a brew before you get started...

I agree that Image::Buffer may be mistaken as a base class for Image:BufferScratch etc. I do like the idea of making it explicit that it is mapped to disk (something that most people are not used to dealing with). If we changed it to something like Image::BufferDisk, is there any reason why we can't include preload as an option in the constructor (either a bool or by passing the strides), and get rid of the BufferPreload class?

Yes, there is a reason: with the standard Buffer, every read/write is done by invoking a function by pointer - i.e. there are two function pointers in the class (one for read, one for write), and these pointers are set to point to the appropriate function in the constructor, based on the type of the data on disk and the type requested as the template argument. These functions will perform the appropriate type conversion, but also any scaling that might have been specified in the header. So there is a definite cost to using them (setting up the function stack, and two floating-point operations). The other two Buffers (BufferScratch & BufferPreload) both access the data inline and natively, with no conversion, scaling or function call - it's a straight pointer dereference, which is as cheap as it's going to get.

What this means, oddly, is that it is more straightforward to merge BufferPreload with BufferScratch... But that makes less sense conceptually than merging it with the standard Buffer. So as things stand, no you couldn't simply have preload be an option on the standard Buffer - the implementations are entirely different.

However, I think there is something we can do about this. Basically, I think we can provide a single Buffer class, providing all the functionality of the above, with a single additional branch point: we just need to check whether data access is direct or via function pointers, and act accordingly (incidentally, the MRtrix 0.2 implementation does precisely this...). This would add a little bit of overhead to data read/write, but probably nothing that can be detected even in the tightest loop. What I'd do is use direct access if the pointer to the data is non-null. Since we're about to use that pointer for the direct operation itself anyway, this can't increase the chances of a cache miss or otherwise slow things down on the CPU (it might for access via function pointers compared to the existing approach, but since it's already slow...). Also, since this test would be inlined, there is a good chance the compiler would optimise it away for the inner-most loop anyway...

With respect to #2, I can see what Rob is getting at. Strictly speaking a buffer is really an image, since it also contains the header info, and it can now be constructed directly from a file name (without creating a header first). The problem is the namespace issue...I agree that Image::Image is a bit weird. Although Image::ImageDisk, Image::ImageScratch, Image::ImageSparse is less weird.

Sorry, I probably should have mentioned that I completely see what Rob was driving at. I actually considered having an Image::Image class myself for the same reasons before we settled on Image::Buffer...

Alternatively we could remove the image namespace (like ITK's image class is just itk::Image). So MR::Image::Buffer would be MR::Image (or MR::ImageDisk, MR::ImageScratch, and MR::ImageSparse), and then all the filters would be MR::Filter::Dilate etc. Since we (currently) don't have filters, interpolation, or registration that work on anything but images (for example meshes), then it's probably not such a big deal removing the Image namespace. Although things like MR::Handler::Base might need the namespace changed to MR::ImageHandler::Base.

I think that's a great idea, and probably the easiest way to sort a lot of this out. You're absolutely right, removing the Image namespace leaves us free to have an Image class, which would admittedly be less confusing than a Buffer. Also, it's pretty clear that all the filters, adapters, etc. would be expected to operate on images unless otherwise stated, so bringing them into the MR::Filter/Adapter/... namespace shouldn't introduce any ambiguity.

When I first started using MRtrix I also found it weird/redundant that I could access the same header info from either a header, buffer or voxel object (although now I just find it useful). At the time, we had to create a header first before creating a buffer, then we would pass voxels around to various functions. Now that we can directly create and pass around buffers (images!), then create voxels whenever needed using auto, does a voxel really need to be able to return buffer/image attributes such as ndim, dim etc?

Short answer: IMO, yes. You don't have to dig too deep in the code to figure out that many of the generic algorithms (not the filters necessarily) operate on voxels, and hence need access to image attributes that would be stored in the header/info. This is particularly important for the whole voxel adapter framework - otherwise you'd have to set up buffer adapters, which I think would be both unwieldy to use and horrible to implement (most of the adaptation relates to or relies on navigating voxels, which is squarely the role of the Voxel class as things stand). And you couldn't set up nested adapters (i.e. an adapter on top of another adapter). It would also make some operations more complicated, for example setting up a simple threaded loop: right now, you can pass a lambda function that just take a voxel or two and operates on that, with access to image attributes if needed. With this alternative framework where voxels don't provide that information, if you did need access to the image attributes, you'd need to set up a proper class constructed from the corresponding buffer/image to glean all the information you'd need to use in the operator - information that would otherwise be readily available as-is.

Also, bear in mind that each voxel needs to store a reference to its corresponding image/buffer to actually access the data; the equivalent methods in the Voxel class are typically thin wrappers around the Buffer/Image version - but crucially not always, notably for adapters. So it makes no sense to actively prevent it from giving you access to that information. What we could do to clarify things is allow access to the image via a dedicated .image() accessor method (or .get_image() if you really like typing), which would make this unambiguous, at the cost of more typing and less genericity (i.e. I could no longer pass a Voxel to a generic algorithm expecting an InfoType or BufferType, I'd need to pass what's returned by its .image() method) - but this would be very difficult to get right for adapters, since the 'image' they operate on is actually a modified version of the original (for instance, the Subset adapter presents a view of a subregion of the image). To make this work, we'd need to have the image() call return a modified copy of the original image, with the right information in it, which would be expensive if some aspect of the relevant information was used in a tight loop. Alternatively, we go for a limited .info() version (as discussed in #210), and each Voxel can keep its own copy of that. But it all seems like a lot of effort to remove something that is currently very simple and convenient, just because it might seem a little unexpected at first.

I stress that while it might be a little bit unexpected, IMO there is nothing ambiguous about it: there's no question about what the information means or which image it might refer to... We could maybe 'fix' this by making the method names a little bit more explicit: it would probably be less jarring if the transform was obtained using a method named image_transform() rather than simply .transform(), and likewise for the .dim(), .ndim(), .vox(), etc. methods. Personally, I think renaming things this way would be a waste of valuable keystrokes for essentially no gain. And if these calls don't match across the Info, Header, Buffer & Voxel classes, they can't all be used interchangeably in a lot of the generic algorithms and functions that we've built up - something that I'd consider a serious regression...

So personally, on that particular front I vote we broadly stick to the status quo.

Possible action plan

That said, having thought about this a bit over the last couple of days, and with the above in mind, I think there are a couple of things we can do. The first few are relatively minor, but if we want to do something serious, I think there is scope for some more major refactoring - as you know, I'm not averse to major overhauls of the code if I think there's a good case for it, it definitely wouldn't be the first time...

Here goes - as usual, comments/suggestions welcome:

remove Image namespace, rename Buffer to Image: this is relatively straightforward, and only involves name changes. It will still be a lot of work, since to do this right we'd also need to move the headers accordingly, and update all the #include lines to match.

merge all Buffer/Image classes into one: harder, but doable. This will change the way the buffer classes are constructed, which may actually be a good thing and sort out another issue that's been bothering me a bit lately: commands often provide the option to include some other image (usually a mask) that can be used in processing. Since there is no way to create an empty buffer, we end up resorting to (smart) pointers to buffers so we can check whether they've been set up later. We also often need to do the same for their corresponding Voxels. This is all a bit messy, it would be simpler to be able to instantiate an empty buffer, with a method to test whether it's usable. This could also be exposed in the Voxel class, so the actual algorithm could just take all the voxels and check whether the optional ones are set or not, rather than take pointers to them to achieve the same effect...

Anyway, what I have in mind is something that would allow syntax like this:

auto vox_in = Image<float>::open (argument[0]).voxel ();
auto vox_in_preloaded = Image<float>::open (argument[1]).voxel (with_strides);
auto vox_out = Image<int>::create (argument[2], vox_in.header()).voxel();
auto vox_scratch = Image<float>::allocate (vox_out, with_strides).voxel();

which can all be done cheaply using C++11 move semantics. Note that in this case, the Voxel class would store a copy of its parent using a shared pointer to ensure it remains in scope, which is probably a better way to handle this anyway. It also simplifies the issues associated with ensuring that the header is not modified after data allocation - more on that later. And it keeps the code really lean, with really no superfluous variables anywhere: just one voxel if that's what you actually need...

Note also that the open(), create() & allocate() calls are static functions, ensuring they can't be used on existing images, which might otherwise be a bit messy to deal with. I'd probably make the default and copy constructors private, and only have the move constructor public.

Allow the creation of empty Image/Buffer: and also allow creation of corresponding voxels. This allows this kind of thing, which at the moment is a bit awkward (using the syntax above):

auto vox_mask = Image<bool>::empty().voxel();
auto opt = App::get_option ("mask");
if (opt.size())
  vox_mask = Image<bool>::open (opt[0][0]);

...

if (vox_mask.valid()) 
  // use it.

Note the use of another static empty() method, specifically for this, which would make sure this is very explicit.

Merge Image::Header and Image::Info: this hasn't been suggested so far, but might remove yet more of the confusion around all of this. We could merge the current Image::Info and Image::Header classes by extending the Image::Info class with a (smart) pointer to an additional structure for all the header-only information (image format, comments, DW scheme, etc). This could be accessed through a header() method, giving access to the information where available (that name might seem odd now, but see the next point). This would allow a single Header/Info class to be used everywhere, since it would be very lightweight (only one additional pointer in it compared to the existing Image::Info class), while retaining the ability to store extended information.

Merge Image::Header/Info into Image class itself: another one that hasn't been suggested before... The main point of keeping the header and buffer separate was to allow loading the header from an image without also loading its data. This is used to provide fast operation for mrinfo for example, but also for other apps where the data of some image might be irrelevant (e.g. mrtransform or tckmap where we specify the template to regrid some other image onto). The idea was that constructing a header loads the header information along with all the information needed to locate the data, and constructing a Buffer would use all that information to actually map / load the data. What we could do instead is trigger a data load only when a Voxel is instantiated. In other words, the Image class is essentially what was the Header before, and requesting a voxel is how we implicitly state our intention to access the data. This is actually a relatively trivial change in terms of the backend, and might make more sense in general. This would mean that:

auto image_in = Image<float>::open (argument[0]); 

would not load the data in any way. image_in can be used where a Header or Info would have been used previously. Then:

auto vox_in = image_in.voxel();

would actually trigger a data load. vox_in would now contain its own copy of the header information, so that image_in can safely be modified or deleted.

Putting it all together

If we go along with all the changes above, we'd end up with just two classes: Image & Voxel (and the image handlers, but they're essentially private anyway). The general idea would be:

What we need to do now

Now if we do decide to go ahead with this, I think it'll need to be done and dusted before release, which doesn't give us a massive amount of time. It's also important we all agree on the details here, this has the potential to get very messy... I'd urge everyone to set up as many tests as is practical for as many parts of the code as they can using the testing framework I added a few days ago, so we can use this to ensure we don't screw everything up in the process...

Let me know what you think, and I'll get started on this if we think it's worth doing once we're all more or less happy with the overall approach...

jdtournier commented 9 years ago

By the way, once you've read the above, there is another minor variation possible, which might make even more sense. Instead of Image & Voxel as the names for the final classes, we could have Header & Image. That means the Header would be closer to what it is now (would only provide access to image attributes, not data), while the Image class would provide all of that information plus the data access methods, as the current Voxel class does. So the Image class really would be a fully-formed image with a bona fide Header, full access to all the data, and copy-constructible for use in multi-threading, etc. Wada y'all reckon...?

draffelt commented 9 years ago

Just started reading... Should have given Piper a bigger bowl of porridge On 14 Apr 2015 2:29 am, "J-Donald Tournier" notifications@github.com wrote:

By the way, once you've read the above, there is another minor variation possible, which might make even more sense. Instead of Image & Voxel as the names for the final classes, we could have Header & Image. That means the Header would be closer to what it is now (would only provide access to image attributes, not data), while the Image class would provide all of that information plus the data access methods, as the current Voxel class does. So the Image class really would be a fully-formed image with a bona fide Header, full access to all the data, and copy-constructible for use in multi-threading, etc. Wada y'all reckon...?

— Reply to this email directly or view it on GitHub https://github.com/jdtournier/mrtrix3/issues/211#issuecomment-92417974.

jdtournier commented 9 years ago

Tell me about it... Took up pretty much the whole day...

On 13 April 2015 at 21:33, Dave notifications@github.com wrote:

Just started reading... Should have given Piper a bigger bowl of porridge On 14 Apr 2015 2:29 am, "J-Donald Tournier" notifications@github.com wrote:

By the way, once you've read the above, there is another minor variation possible, which might make even more sense. Instead of Image & Voxel as the names for the final classes, we could have Header & Image. That means the Header would be closer to what it is now (would only provide access to image attributes, not data), while the Image class would provide all of that information plus the data access methods, as the current Voxel class does. So the Image class really would be a fully-formed image with a bona fide Header, full access to all the data, and copy-constructible for use in multi-threading, etc. Wada y'all reckon...?

— Reply to this email directly or view it on GitHub <https://github.com/jdtournier/mrtrix3/issues/211#issuecomment-92417974 .

— Reply to this email directly or view it on GitHub https://github.com/jdtournier/mrtrix3/issues/211#issuecomment-92490131.

Dr J-Donald Tournier (PhD)

_Senior Lecturer, _Biomedical Engineering

Division of Imaging Sciences & Biomedical EngineeringKing's College London

A: Department of Perinatal Imaging & Health, 1st Floor South Wing, St Thomas' Hospital, London. SE1 7EH T: +44 (0)20 7188 7118 ext 53613 W: http://www.kcl.ac.uk/medicine/research/divisions/imaging/departments/biomedengineering http://www.kcl.ac.uk/medicine/research/divisions/imaging/departments/biomedengineering

draffelt commented 9 years ago

Good ideas on all fronts! I agree with all of them. This seems much cleaner and easy to understand from a new user perspective. The idea to explicitly create an empty image is a good one. Mask ptrs are always a little messy.

I also like the idea of renaming Image & Voxel to Header & Image. Especially since you already mentioned above that an Image (ie a Header) class "contains all the basic information to represent an image". It also solves the issue of a Voxel having image attributes ;)... not that I really had any issue with this.

So just to confirm, we can load and access an image like this:

auto image_in = Header<float>::open (argument[0]).image ();
for (auto i = LoopInOrder() (image_in); i ++i)
    image_in.value() *= 2;

I will start to write some test apps over the next few days.

jdtournier commented 9 years ago

Cool, glad you agree. I'll get cracking with this asap.

By the way, the example you've got there is almost right. It would compile, but crash at runtime since you're trying to write to a read-only (input) image... But in terms of syntax, yes that would be the idea.

jdtournier commented 9 years ago

OK, got started with this. It's incredible how much this seems to simplify the header. One of the immediate benefits is that all member variables can now be made public, since they are modifiable by default, and the Voxel / Image class would store a const copy. There's a lot of cruft that can immediately go, in terms of all the accessor methods, etc - not to mention this ConstHeader class, which was designed as a semi-modifiable version for use in the Buffer classes...

However, one question it does throw up is whether we need to stick to the dim(), ndim(), vox() & stride() methods. What these do is simply act as thin wrappers around a std::vector of Axis classes, which are themselves really simply objects with dim, vox & stride members. We could omit these accessor methods, which would mean changing from this:

H.set_ndim (3);
H.dim(0) = 10;
H.vox(1) = 1.2;
H.stride(2) = -1;

to this:

H.axes.resize(3);
H.axes[0].dim = 10;
H.axes[1].vox = 1.2;
H.axes[2].stride = -1;

So this would entail a bit more typing, but I think it might actually be more legible this way. What do people feel about this? Happy to keep the old accessor methods if it's felt to be superior, but I'm actually leaning towards the slightly more verbose second option right now...

Lestropie commented 9 years ago

This is all a bit much for my walking-wounded brain right now, but I'll add what I can.

Ah crap. This has now turned into a full blown mutiny...

Don't worry, if there's a mutiny you'll know about it... you won't be able to log in to GitHub :-P

Basically, I think we can provide a single Buffer class, providing all the functionality of the above, with a single additional branch point: we just need to check whether data access is direct or via function pointers, and act accordingly (incidentally, the MRtrix 0.2 implementation does precisely this...). This would add a little bit of overhead to data read/write, but probably nothing that can be detected even in the tightest loop. What I'd do is use direct access if the pointer to the data is non-null. Since we're about to use that pointer for the direct operation itself anyway, this can't increase the chances of a cache miss or otherwise slow things down on the CPU (it might for access via function pointers compared to the existing approach, but since it's already slow...)

My concern here would be the breaking of the on-die pipelining due to the additional branch. But won't know unless we have the code to benchmark.

remove Image namespace

Although I agree that there's a lot of the stuff in MR::Image that can't be applied to anything that's not an image, this change would mean that there'd be a lot of stuff in the MR namespace that operates on images and a lot of other stuff that doesn't.

Anyway, what I have in mind is something that would allow syntax like this:

Another worth considering here is the BufferAtomic: it's sitting there in a branch ready to merge, it just doesn't have an application yet. This would either need to use get/set functions with atomic casts (and hence use the get/set functions even though the scratch data is set), or allocate an array of atomics; either way, more branching.

Also: Having Buffer / BufferPreload / BufferScratch means that BufferSparse sits relatively sensibly in the scheme of things, given how different the storage mechanism is. If the first three are homogenised, the sparse stuff may be harder to find a neat place for.

However, one question it does throw up is whether we need to stick to the dim(), ndim(), vox() & stride() methods.

I've found a couple of places where calling set_ndim(3) resulted in undefined behaviour, and sanitise() needed to be called; keeping these as functions may give a mechanism to detect when such sanitisation is necessary...? But I think you're right in that functionalising these was to force const-ness in certain instances, which may no longer be required if the accessor class stores a const copy.

draffelt commented 9 years ago

Personally I prefer the idea of keeping the accessors methods and hiding the details of how these attributes are stored in a vector of axes. I'm not sure why having to type axes makes it more legible. I agree that these attributed are all related, however it feels like you are just making them public because you can, without any gain from the users point of view. For example it is pretty obvious what set_ndim does. And I like to type as little as possible.

In any case, I don't feel that strongly about it.

jdtournier commented 9 years ago

OK, before I start responding to some of the comments, it's just occurred to me that the syntax I'd suggested is probably not quite right. I think the Header should not be a template class, only the Image. So to take @draffelt 's example above, the syntax will probably end up like this:

auto image_in = Header::open (argument[0]).image<float> ();

rather than what I'd suggested before:

auto image_in = Header<float>::open (argument[0]).image ();

i.e. it's the .image() call that will be templated, not the Header returned by open()/create()/allocate(). This keeps things closer to the way they're handled now.

OK, now onto the other comments:

@Lestropie :

My concern here would be the breaking of the on-die pipelining due to the additional branch. But won't know unless we have the code to benchmark.

Pretty sure that wouldn't be a problem - modern CPUs are pretty aggressive with their optimisations. In any case, I reckon the compiler itself would optimise some of that away... We can probably help by making sure the condition variable is declared const (the pointer's address, that is, not the pointed-to data) - hopefully that'll give it a hint that it's not going to change from one iteration to the next, so the branch only needs to happen once before the loop even begins...

Although I agree that there's a lot of the stuff in MR::Image that can't be applied to anything that's not an image, this change would mean that there'd be a lot of stuff in the MR namespace that operates on images and a lot of other stuff that doesn't.

I'm not sure there's all that much, really. Aside from some of the string helper functions, exceptions, progressbar, and the timer, there's not a lot else (?). Just about everything should be in its own namespace.

Also: Having Buffer / BufferPreload / BufferScratch means that BufferSparse sits relatively sensibly in the scheme of things, given how different the storage mechanism is. If the first three are homogenised, the sparse stuff may be harder to find a neat place for.

Yes, this could be an issue. But thinking about it, all it would take is to add the equivalent sparse static calls to get the Image (i.e. what was the Voxel) - so in addition to Header::open()/create()/allocate(), you'd also have Header::open_sparse()/create_sparse()/allocate_sparse(). These would return a SparseImage object instead of the normal Image. Within that object, you'd be free to add a shared_ptr to any additional data that might be needed for sparse image access. I reckon it would actually be relatively straightforward.

Another worth considering here is the BufferAtomic: it's sitting there in a branch ready to merge, it just doesn't have an application yet. This would either need to use get/set functions with atomic casts (and hence use the get/set functions even though the scratch data is set), or allocate an array of atomics; either way, more branching.

I think we could adopt the same approach as above. Or better yet, use template specialisation or SFINAE to return a modified Image class that implements the atomic operations when the value_type is std::atomic<>. Either way, it can be done and can probably re-use a lot of the existing code.

I've found a couple of places where calling set_ndim(3) resulted in undefined behaviour, and sanitise() needed to be called; keeping these as functions may give a mechanism to detect when such sanitisation is necessary...?

Well, really sanitise() should be called within open()/create()/allocate() to finalise the header - the header shouldn't be assumed valid before then, and shouldn't be modified / modifiable after that. I was actually thinking of making it private, although I can see no reason to prevent its use. If you get undefined behaviour, it would indicate some other invalid assumption somewhere. Maybe you can give me a specific example if you think this might indeed be a problem.

But I think you're right in that functionalising these was to force const-ness in certain instances, which may no longer be required if the accessor class stores a const copy.

Exactly. Really, the way I see it, the Header class is now essentially a glorified C struct, which you can manipulate any way you want - it has no side effects. Once you do call open()/create()/allocate() though, it shouldn't be modified. I'll probably enforce this by getting these functions to return a const Header. The empty() call can return a non-const version, and the copy constructor can be used to copy an existing const header and modify it.

@draffelt :

Personally I prefer the idea of keeping the accessors methods and hiding the details of how these attributes are stored in a vector of axes.

Yes, I'm in two minds about that one. I like the idea of keeping this class as close to a bare-bones C struct as possible, to convey the sense that it is really just a placeholder for information, which will be fully sanitised and checked when the data are being accessed.

But the reality is that removing these methods would require extensive changes throughout the codebase to match the new syntax... Let's keep them in.

jdtournier commented 9 years ago

Another minor issue these changes brings up:

Currently, you would use the Voxel class's operator[] methods to get/set the voxel's position along each axis. This kinda made sense when the class was called Voxel - but if we now essentially rename it to Image, that would probably seem pretty weird...

I think we need to change to syntax for that to something a little more explicit, so for example:

vox[0] = 0; 
++vox[1];
vox[2] += 10;

would change to:

image.pos(0) = 0;
++image.pos(1);
image.pos(2) += 10;

Any objections to this...?

Lestropie commented 9 years ago

Hmm. Might feel weird copy-constructing an 'image' for multi-threading...?

jdtournier commented 9 years ago

Yes, that had also occurred to me... But in actual practice, I don't think you'd notice much, if at all. You won't generally be copy-constructing explicitly, you'd have an Image as a member of some Functor class that you'd just pass to the appropriate Thread call, and that call would create copies behind the scenes, so I don't think that would feel weird at all. Within the Functor itself, you'd likewise just be operating on the Image, and that would probably feel quite natural. So I reckon that won't actually be a problem in practice. What do you reckon...?

Lestropie commented 9 years ago

You won't generally be copy-constructing explicitly, you'd have an Image as a member of some Functor class that you'd just pass to the appropriate Thread call, and that call would create copies behind the scenes, so I don't think that would feel weird at all.

I tend to be explicit in my copy-construction; probably too many bad experiences. But would probably get used to it.

sanitise() ... Maybe you can give me a specific example if you think this might indeed be a problem.

I think tckgen was one recently... If you used a 4D image as the template image, tckgen would reduce it to 3D, but it'd give some pretty unusual behaviour down the track (was a bugger to debug from memory).

My concern here would be the breaking of the on-die pipelining due to the additional branch. But won't know unless we have the code to benchmark.

Pretty sure that wouldn't be a problem - modern CPUs are pretty aggressive with their optimisations. In any case, I reckon the compiler itself would optimise some of that away... We can probably help by making sure the condition variable is declared const (the pointer's address, that is, not the pointed-to data) - hopefully that'll give it a hint that it's not going to change from one iteration to the next, so the branch only needs to happen once before the loop even begins...

I wonder if it would be possible to use a latent template variable to SFINAE out the function call? Would be a bastard without auto / decltype, but with extensive use of those...?

jdtournier commented 9 years ago

I tend to be explicit in my copy-construction; probably too many bad experiences. But would probably get used to it.

Nowadays, I prefer to explicitly ask for the default copy constructor - unless the default copy constructor of one of the members genuinely behaves differently from what I need, and its behaviour is sensible in other contexts, so it isn't just a matter of fixing that member's copy constructor. I find this to be much easier to follow when looking at the code, I don't need to figure out whether there's a bug or an omission in the copy constructor... I find I very rarely need to implement my own copy constructor. A lot of this relies on proper use of the smart pointers - this is why we have RefPtr (soon std::unique_ptr) and Ptr (soon copy_ptr), which covers most use cases for this.

I think tckgen was one recently... If you used a 4D image as the template image, tckgen would reduce it to 3D, but it'd give some pretty unusual behaviour down the track (was a bugger to debug from memory).

Do you mean tckmap...? In any case, I'm not sure I see why there would be a problem here. Assuming you don't have a working example handy, let's worry about this nearer merge time, we should be able to test for this in the regression testing we'll need to do anyway.

I wonder if it would be possible to use a latent template variable to SFINAE out the function call? Would be a bastard without auto / decltype, but with extensive use of those...?

Yes, I was thinking this might be possible. All it would need is for the get_image() calls to return an Image with the template parameter set one way or the other, depending on whether the Preload option is given or not (which would be achieved using plain old function overloading). But the problem is I'm not sure it would necessarily be desirable. To expand: Using the single Voxel/Image class, and branching on the presence of a direct access pointer, we have the opportunity to optimise away the type conversion function calls for non-preloaded, native data type, unscaled images - which will actually be a very common occurrence. So if you pass a float32 image with scale 1 and offset 0, the backend would be able to memory-map that and access the data directly using this framework, which may speed up a number of applications where preloading isn't necessarily useful (i.e. most single-pass applications). SFINAE relies on compile-time tests, so you couldn't do this run-time optimisation - it would either be preloaded and direct-access, or not preloaded and function call access (which is the current situation anyway). I'm thinking the opportunity for run-time optimisation might be worth capitalising on - assuming the cost of the branch is negligible...

draffelt commented 9 years ago

I don't think using image[0] = 10 is that weird to change the index.

I'm not that opposed to image.pos(0) either, however for some reason the word position makes me think of a continuous variable (ie scanner space position) and not a discrete voxel index.

What about if we could have something like image.vox[0] = 10. We could then access the voxel size using image.vox.size[0]. image.vox(0)has always struct me as not immediately obvious it returns the voxel size.

While I'm at it ;) As we have discussed before, I think image.dim(0) is also less obvious than image.size(0). I know you could say to someone "what are the dimensions of the image?". However it's probably more common to say "what is the image size?". Also we have mrresize as a command. Gimp also asks you for the 'image size' when recreating a new image. ITK also uses the word size http://www.itk.org/Doxygen/html/classitk_1_1Size.html

jdtournier commented 9 years ago

OK, happy to go with image.vox[0] = 10 syntax. Also happy to use image.size[0] for what was dim(0) - I agree it probably makes more sense. Not so sure about the replacement for vox(0) - I agree it wasn't not all that clear, but it was nice and short... Going for image.vox.size[0] would be really weird to implement, but maybe image.voxsize[0]? Also wondering about ndim() - never really liked it all that much, but I just can't think of a better term that is short and unambiguous...

jdtournier commented 9 years ago

Actually, thinking about it a little bit more, if we're happy to stick with image[0] = 10 for navigation, let's do that. Otherwise it'll create total havoc as we try to refactor the code, and forget whether that vox should be replaced by voxsize, or has it just been added in...

draffelt commented 9 years ago

Cool. Happy with that. I've also got no objections against image.ndim()

jdtournier commented 9 years ago

Sorry - more thoughts...

How about image.index(0), rather than image.vox[0]? It's probably clearer still, and avoids confusion in the transition period. Also, I'd avoid square brackets, it'd make things too messy implementation-wise...

So we'd end up with something like:

draffelt commented 9 years ago

Sounds good. I was going to suggest index. I think it makes it very obvious what is being modified.

jdtournier commented 9 years ago

Sounds good. I'll wait before going ahead in case anyone else has a strong opinion on the matter.

jdtournier commented 9 years ago

OK, further comments on the coming changes:

I'll push a branch with the current state of things when I have a minute. What I'm worried about most at the moment is merging back to master if significant changes have taken place in the meantime... I'll merge most of the pull requests we have pending now to make my life easier later - most of the pull request seem relatively straightforward (aside from the ODF display - working on that now).

draffelt commented 9 years ago

All sounds good to me. Looking forward to testing it out. On 21 Apr 2015 7:19 pm, "J-Donald Tournier" notifications@github.com wrote:

OK, further comments on the coming changes:

-

we'll stick to using access functions rather than accessing the members of the Header directly - this allows derived classes, such as the Image and any Adapters, to provide an identical interface to the Header and to each other, which is a Good Thing for use in generic algorithms.

I'm thinking of removing from the Header class anything that isn't strictly necessary for operation. This includes comments and the DW scheme, both of which can easily be stored as generic entries in the Header anyway. Currently, generic entries are stored in the std::mapstd::string,std::string that Image::Header derives from. The changes I'm planning would stored this in the same structure, but it would be a member rather than a base class, accessible via the .keyval() method. Getting at the comments and DW scheme would be a matter of parsing the relevant text entries, and shouldn't require a great deal of modification to the existing code since we should be using DWI::get_DW_scheme() or DWI::get_valid_DW_scheme() anyway...

there's a whole bunch of other implications that go along with the changes proposed, but they're relatively minor and too many to list here. Hopefully they'll make sense once they get merged...

I'll push a branch with the current state of things when I have a minute. What I'm worried about most at the moment is merging back to master if significant changes have taken place in the meantime... I'll merge most of the pull requests we have pending now to make my life easier later - most of the pull request seem relatively straightforward (aside from the ODF display - working on that now).

— Reply to this email directly or view it on GitHub https://github.com/jdtournier/mrtrix3/issues/211#issuecomment-94717467.

jdtournier commented 9 years ago

Closing, and moving conversation onto pull request #225.