Scikit-image: Proposal: deprecate `multichannel=` in favour of `channel_axis=`

[jni]

This might need a SKIP. =P But won't if we get unanimous agreement!

The name multichannel is a bit clunky. Worse than that, it does not quite provide sufficient information: we hardcode the channel axis as -1, which might actually be incorrect. Allowing for arbitrary channel axes allows (1) giving the data in acquisition axis order, which often does not have the channels as the last axis, and (2) allowing channels to be the leading axis, which in many cases is more efficient, as the processing happens independently on each channel (e.g. gaussian filtering), and this is more cache friendly if each channel is contiguous.

Over at napari, we've settled on channel_axis=None being equivalent to multichannel=False, and channel_axis=<int> to point to an axis to interpret as channel rather than spatial data. I quite like it, and seeing an issue like #4293, where we might want to either specify no channels or a different channel axis makes me think that we can refine scikit-image to use this also.

Looking forward to hearing your thoughts!

[grlee77]

I also happened to use exactly this channel_axis naming in some local code dealing with multichannel MRI data, as some software assumed channels last while others assume channels are first.

I am +1 on moving toward using channel_axis, but how do you anticipate handling existing code that assumes channels are the last dimension? The simplest solution is probably to manually call swapaxes on function input and then undo the swap on output, but that comes with some performance penalty for cases where the channel axis is not already last.

[jni]

how do you anticipate handling existing code that assumes channels are the last dimension? The simplest solution is probably to manually call swapaxes on function input and then undo the swap on output

Yes, that would be my proposal. The API matters much, much more than performance. We can deal with performance later, when it becomes a measured problem. This would not be a performance regression for any existing code. We can do plenty of approximate tweaks, like, if channel_axis is 0, and the image size is greater than 4MB (or your L1 cache, or whatever), then do a copy of the axis swapped-data before processing. But we should not clutter the code with these until someone complains.

[lagru]

Not sure about the feasibility of the optimizations but I really like API a channel_axis argument would enable. It kinda reflects NumPy's well established axis argument. +1

[hmaarrfk]

My biggest fear is that the release cycle is too slow to make all this maintenance work worthwhile.

It seemed that the biggest fear from our last conversation was that the deprecation cycle would be thrown off if we were to move to more frequent releases. That said, slow releases mean that we don't get feedback about new features for a long time.

These big swooping changes require lots of input from core devs, which meqn that we often don't have time to evaluate other PRs that work on content and new analysis techniques.

TLDR, I am +1, but would be +2 if combined with a more frequent release plan.

[emmanuelle]

Also +1 on this.

[hmaarrfk]

What about cases where there are multiple channels??? Fluorescence and z for example?

[jni]

z is not a channel, it’s a spatial axis.

[soupault]

@hmaarrfk I was thinking the same. :) We should probably better define what image is and the extents of our scope in relation to such definition. For example, do we consider 2D Fourier transform (two frequency axes) to be a valid image for the library (i.e. support it as a valid input for most of the functions)? What about the data array where we have two dims of spatial data, third axis is frequency, fourth is phase, fifth is time? Should we support such data, in general? If so, we may want to consider adding phase_axis= and time_axis=. That would make us closer to xarray, since the spacing for the special dimensions is often not constant and would need axis indices.

[hmaarrfk]

Fourier transforms: Are complex, much of the C code doesn't work with complex numbers, and cython doesn't support them well, at least the last time I tried it didn't.

What kind of phase are you speaking of? Is phase not a channel?

I wouldn't worry too much about what the channels represent, more that there are "dimensions we want to operate on, and dimensions we want the operation to be repeated over".

[soupault]

What kind of phase are you speaking of? Is phase not a channel?

Hmm... I typically associate channel with frequency or frequency range. Although, I guess, you are right, phase is a channel, but can be a _second_ channel axis.

[hmaarrfk]

So you would have amplitude and phase for each wavelength of light you are measuring?

[jni]

@hmaarrfk @soupault I think you are applying a lot of scope creep to this issue. I propose the following, and only the following:

---

• multichannel=True becomes channel_axis=-1.
• multichannel=False becomes channel_axis=None.

The proposal allows us to maintain backwards compatibility with all our current features and express one new thing, which is the case where the channels are arranged along a different axis, say, the leading axis, with:

• channel_axis=0,
• or, for say TZCYX volumes (common in microscopy), channel_axis=2.

This is straightforward to implement (just needs a NumPy transpose before and after existing code), and is strictly more expressive than what we have now, and imho results in a more readable API.

---

Metaphysical discussions about what constitutes a channel are interesting but imho only very mildly affect the discussion. I just met someone (@ericpre) who had, at every pixel of the image (spatial location in a sample), a full X-ray diffraction pattern (ie another 2D image). In this case, you could either ravel the pattern, or you could have a tuple input to channel_axis=, but that can come later, the same way that axis= gained the ability to receive tuples in NumPy ufuncs.

Typically, pixels have spatiotemporal extent, while channels are measurements at the same location. When displaying images, channels are usually overlaid on top of each other. Other axes, like time, z, etc are sliced. Sure, there are contexts when you might want to treat time as a channel, but they are in the minority, and again, you can use channel_axis=a_tuple for those situations. To go back to TZCYX, you would pass in channel_axis=(0, 2), e.g. if you wanted to apply a Gaussian filter over the space dimensions. (You could also do channel_axis=2, sigma=(0, 3, 5, 5) to prevent smoothing over the time axis.)

And I'm not sure why complex dtype images would affect the discussion at all. It's an orthogonal issue.

[soupault]

@jni My point above is rather an attempt to understand whether we should use this opportunity and make things a bit differently / more flexible. I generally agree with your definition of channels. However, I believe there are cases with multiple channels where one does not want to consider them similarly during the processing (e.g. filtering, inpainting). One of the examples could be a dual-energy X-ray captured in the form of both the absorption and the phase contrast images. On the other hand, such cases with multiple channels are probably pretty rare and are easily solvable by just iterating over the proper axis. To assist this, we should probably document how the channels are handled in each of the functions with multichannel support.

Sorry, I did not make that clear that I support your original proposal. I think having a tuple as a valid parameter value in channel_axis would also be great. And I did not mean to block the proposal in any way.

[grlee77]

multichannel=True becomes channel_axis=-1.
multichannel=False becomes channel_axis=None.

I am +1 on this and also on potential extending channel_axis to a tuple in the future as well.

An example where a tuple of axes would be easy to support is in denoise_wavelet, as recent releases of PyWavelets accept a tuple of axes to perform the transform over (in this case the transform axes would be all of those NOT in channel_axis).

[hmaarrfk]

Jni, happy to keep the discussion contained.

I see two things being good to have as part of this change:

1. The introduction of the channel_axis parameter on an as needed basis.
2. The inclusion of a note on which axis is likely to have the best performance.

I'm +1 on this

[jni]

@hmaarrfk what do you mean "as-needed"? Do you mean a slow migration function by function? Or...? I would prefer a single one-shot PR, as I think it would be confusing and inconsistent to have both styles in the library...

[hmaarrfk]

I think it will take alot of effort to do this right.

Tests will have to be written, and reviewed. Doing it all at once seems impossibly hard

[jni]

@hmaarrfk challenge accepted. =D

[rfezzani]

@jni, this looks like a Meta PR! With the functions cross calls, the tests... I agree with @hmaarrfk, it will be hard to do end even more harder to review...
I know that I am may be not the right person to say that (I wrote some huge PR :stuck_out_tongue_closed_eyes: ), but this one will be gigantic!
Why not migrate module by module, if function by function is too slow?

[grlee77]

Why not migrate module by module, if function by function is too slow?

Making PRs on a module-by-module basis sounds reasonable to me. I would be willing to help cover a couple of modules if we split the effort.

I don't think it makes sense to have a new release with only a partial conversion, but in between releases, I don't see how it hurts to convert module by module in the development branch.

[rfezzani]

I just realized that the deprecation of multichannel in favor of channel_axis can be easily made using the deprecate_kwarg decorator that I define in #4158 (still needs to be merged :wink:).

[jni]

@rfezzani I thought of that, but it's not quite so simple: channel_axis=0 cannot be expressed currently, and requires a transpose of the input and a transpose back of the result (at a minimum — in some cases, that will not be the most efficient way to do it).

[rfezzani]

This decorator can be adapted to manage this: it currently only defines kwarg_mapping for describing the mapping between the old and new argument names, but we can add the kwarg_values to define the default value of the new argument name. What do you think about this solution?

[hmaarrfk]

I find decorators to be cool software tricks, but difficult to grasp if you aren't software savvy

[rfezzani]

The proposed decorator is for internal use only. I already listed its advantages here https://github.com/scikit-image/scikit-image/pull/4250#issuecomment-544980145.
Even if in #4250 the

Finally, it is really easy for contributor's to understand what is going on without decorators.

argument trump everything, to me,

less noise introduced to the documentation, since the deprecated argument is replaced by the new one instead of adding a new argument to the argument list, and only the users concerned by the deprecation are impacted.

argument is a major one...

[hmaarrfk]

Is the decorator meant to be temporary or permanent?

[hmaarrfk]

I would really want https://github.com/scikit-image/scikit-image/issues/4454 to be accepted before moving on this deprecation.

I think that this deprecation is going to be insanely invasive and cause many headaches. I think it is important for us to understand the magnitude of the change, before we push it on to others.

[rfezzani]

The decorator is made to be part of the deprecation cycle. The deprecation cycle starts when a function is decorated and ends when the decorator is removed. So yes, clearly, the decorator is temporary.
I am :+1: for #4454

[sciunto]

xref: #2613 default values of multichannel argument

[grlee77]

xref: #2613 default values of multichannel argument

Thanks @soupault. That list was getting pretty outdated, so I just updated it (at the top of the issue)