PyDIP, where should it go from here? (pt 2)

[crisluengo]

Issue #25 was getting too long, I added a recap at the end. Here's a clean slate to continue the process.

• We need lots of usage, lots of feedback from users, and lots of time polishing the interface.

• But ideally we'd take this to a different level, in a way similar to what DIPimage is. For example, we could write a GUI as exists in DIPimage. And maybe make the module more "Pythonic" by changing function names.

• We also need to find a way to automatically generate documentation from Doxygen. One approach is here.

[crisluengo]

dip.Polygon is rather useless. It has some methods to measure the polygon, and it can be converted to a NumPy array to get its vertices. It'd be nice if it had __getitem__, __len__, __iter__ and so on, so it can be used as a list. (this is complete)

[grlee77]

One issue I found confusing when trying out PyDIP, was that I had to reverse the order of elements in tuples provided as arguments from what I would have naively expected. I give two concrete examples below:

example 1

Assume I want to do Gaussian filtering of a 3D volume via dip.Gauss(image, sigmas) where sigmas is a 3-tuple of floats. I would expect that specifying sigmas = (1.0, 2.0, .4.0) would smooth the first axis of numpy.ndarray image with sigma 1.0 and the last axis with 4.0, but it is actually the opposite which happens. To get the expected behavior I have to specify dip.Gauss(image, sigmas[::-1]). I assume this is related to the statement near the bottom of this page:

The mapping from images to arrays causes the indexes to be reversed: The first array index corresponds to the last image index, and vice versa. If an image is indexed as img[x,y,z], the corresponding array is indexed as array[z,y,x].

example 2

Suppose I want to perform morphological operations with a rectangular structuring element of shape (5, 9). In scikit-image, I could define the structuring element via footprint=np.ones((5, 9)) and call skimage.morphology.binary_erosion(image, footprint).

However, in DIPlib, to get identical output with a rectangular structuring element I would need to create the structuring element using dip.SE((9, 5), 'rectangular') (note that the shape must be specified in reverse order!).

proposed solution

Should these function's wrappers detect if the input is a NumPy array and reverse the entries of tuples like this automatically? That would seem to be a more user-friendly approach. Otherwise, I suspect it is likely that many users of the library will make these types of mistakes. I hope this feedback is helpful to the team!

[grlee77]

In case it is of interest, I came across the examples above when I was prototyping possibilities for offloading scikit-image computations to different backends. There are some cases where the behavior of equivalent routines across scipy.ndimage and DIPlib is identical to within floating point precision (e.g. Gaussian filtering, morphology). In those cases, we could optionally dispatch to DIPlib instead of scipy.ndimage when available and get a performance benefit.

There are some examples and benchmarks for a small number of functions in this demo project's README. You will see that DIPlib beats vanilla scipy.ndimage or use of parallel blockwise processing via Dask+scipy.ndimage, losing out only to a couple of GPU-based cases.

[crisluengo]

@grlee77 Thanks for your thoughts! I see your point, and agree it would be convenient if things matched up.

I'm not sure if we can make your proposed solution work without significant changes to the bindings. We currently take as input any object with a buffer protocol, the bindings automatically create a dip::Image object around it, and pass that to the C++ function. The C++ function currently doesn't know that the input wasn't a dip.Image object. If we are able to make that distinction, we would have to add a bit of code to each function to reverse the relevant array parameters. Currently the bindings for most functions just indicate to the compiler which C++ function to make available in Python, most of the actual code is for translating Python types to DIPlib types. So adding a bit of code to each function would probably mean multiplying the size of the code for the bindings by 5 or 10 or so. This would be a significant amount of work and make it a lot harder to maintain the bindings and verify that they're correct. It's not something I'm keen on.

Though there might be some simple solution that I can't think of right now. Will keep thinking about this!

Here's my reasons for reversing the dimensions of the NumPy arrays:
In image processing, dimension order has always been (x,y,z), and so that is also the ordering that DIPlib has used since before it was popular to do image processing in MATLAB or Python. But NumPy indexes the array as (z,y,x), and that order is respected by skimage and imageio and PIL and matplotlib and so on. And so the DIPlib Python bindings had to take arrays representing an image in (z,y,x) order and make it useful in a C++ library that uses normal dimension ordering. If we don't reverse the dimensions of the array, images loaded through imageio, PIL, OpenCV, ... will all be seen sideways in DIPlib's viewer, and images loaded through DIPlib will be seen sideways in matplotlib, and saved sideways to file.

[wcaarls]

It seems the main problem is import/export and visualization. Can we make a global option that instructs those functions to interpret the axes differently? For the viewer at least it’s just a matter of changing the default axes and their names.

Of course, the fact that numpy arrays would have non-normal strides when converted to images and vice-versa would probably affect performance.

[crisluengo]

There are some DIPlib functions that process the image in order of coordinates (i.e. increasing x first), but most functions should work equally fast no matter what the dimension order is. Remember that in MATLAB we are always dealing with images with non-standard strides, so this has always been important to get right.

So this would require:

1. BufferToImage and ImageToBuffer to reverse sizes and strides arrays optionally, depending on a global variable.
2. Add a function that can be used to change that global variable.
3. Change PyDIPjavaio.ImageRead*, PyDIP_bin.ImageRead* and PyDIP_bin.ImageWrite* to reorder dimensions after calling the C++ function (or before in the case of writing) depending on the global variable.
4. Change PyDIPviewer.Show, and some methods and properties of PyDIPviewer.SliceViewer, to reorder dimensions depending on the global variable.

Am I missing something?

This seems like a very doable proposition.

[wcaarls]

Good to know about the performance! I was hoping that would be the case.

Do you want the option to be specific to PyDIP, or also available for users of the C++ library? In the latter case, the C++ functions themselves would do the translation.

[grlee77]

Thanks for looking into the possibilities here. It sounds like perhaps there is a solution that won't be too difficult from an implementation and maintenance standpoint. If it turns out to be harder than expected, just making an explanation of the current behavior more explicit in the PyDIP docs would help.

[crisluengo]

@wcaarls I was thinking of having this specifically in the Python interface. Didn't think about C++. Is there a use case for it in C++? Also, that would require adding a global variable to DIPlib, which I'm not keen on (I removed all the globals that DIPlib 2 had, we only have the OpenMP number of threads as a per-thread global variable right now).

[wcaarls]

@crisluengo No use case that I know of, but it impacts how I implement it in the viewer. If it is PyDIP only it needs to be a viewer option, whereas if it is a DIPlib variable the viewer can query it by itself. I'll make a viewer option then.

[crisluengo]

@wcaarls It should be easy to reorder dimensions in the function in the bindings, see dip::Image::PermuteDimensions. This is a cheap operation that doesn't copy data. You shouldn't need to touch anything outside PyDIP to do this.

dip::Image tmp = input; // because input is likely a const&, we can't and don't want to change it
dip::uint ndims = tmp.Dimensionality();
dip::UnsignedArray order(ndims);
std::generate(order.begin(), order.end(), [& ndims]{ return --ndims; });
tmp.PermuteDimensions(order);

Maybe we can add this as a member function for dip::Image, to make it a simple call. Something like dip::Image::ReverseDimensions.

[crisluengo]

@grlee77 I've added some more explicit description of what happens when you mix NumPy with DIPlib. Do you think it's clear enough now? b07d898

[grlee77]

@grlee77 I've added some more explicit description of what happens when you mix NumPy with DIPlib. Do you think it's clear enough now? b07d898

Yeah, I think that helps. (scikit-image should be lowercase, though)

[crisluengo]

@wcaarls I have pushed the changes for dip.ReverseDimensions(). I have left most of our alternative solutions in this branch for future reference: https://github.com/DIPlib/diplib/tree/optional_reverse_dimensions (I think only your cached module version is not there, feel free to add it).

@grlee77 In the next release, calling dip.ReverseDimensions() will switch the dimension order to match numpy and skimage. This function is explicitly designed to not allow you to go back and forth, and to avoid weird results you should call it immediately after loading the diplib module.

[crisluengo]

I'd like to change the * and *= operators to use dip::MultiplySampleWise() instead of the dip::Multiply(), and then implement __matmul__ and __imatmul__ to use dip::Multiply(). But I'm afraid of introducing such a breaking change.

I guess we can tell people to do

dip.Image.__mul__ = dip.Image.__matmul__
dip.Image.__imul__ = dip.Image.__imatmul__

to revert to the old behavior? Like I did with the dip.Histogram change.

The Python bindings have lots of awkward things that are hard to fix without breaking existing scripts...

@wcaarls What do you think?

[wcaarls]

@crisluengo Although inconvenient, I think pydip should be treated as a moving target for now, since we don't have a clear idea of where it should go. So I don't think we should get hung up too much on breaking changes, as long as they are a definite improvement and are announced clearly. I think the proposed change qualifies.

[crisluengo]

@wcaarls I think that's a good policy, it will make it easier to make improvements.

I've been doing version numbers like this: 3.x.y

• If y increases it means you can replace libDIP.so without recompiling the program that uses it.
• If x increases it means you must recompile the program that uses it, but you don't need to make changes to the code.

Obviously this refers to the C++ library only, not anything else.

Would it be OK to increase only y if we make a non-compatible change to PyDIP like the above? Or would it be doubly surprising that going from 3.2.0 to 3.2.1 your Python script no longer works as previously?

It's kinda awkward having a single version number for all these different components. I thought it was convenient putting it all in the same project, but maybe it's not?

[wcaarls]

Do you see the 3.x.y as major.minor.patch (semver.org) or epoch.major.minor? In general, people will assume the first.

Theoretically, breaking changes require a major version bump, but we can treat PyDIP as a de facto development release and only bump the minor version. I do think the patch version (if any) should be reserved for bugfixes.

Separate component versions will probably be more confusing. I don't think 3.42.0 is too problematic.

[crisluengo]

I've been doing {incompatible}.{API compatible},{ABI compatible}, but that is explicitly related to the DIPlib binary. And we haven't really had a 3.x.1 yet, so I don't think I got to explain the version scheme anywhere yet.

I do think the patch version (if any) should be reserved for bugfixes.
I don't think 3.42.0 is too problematic.

Let's do it that way then. Whatever is least surprising. :)

I'll add some words on the website and README to indicate PyDIP is still in development and might change in incompatible ways. People should pin their version number and read release notes before upgrading.

Thanks!

[crisluengo]

I've been pointed to the __array_wrap__ method, used by NumPy when doing operations on foreign objects, so it can return an object of the same type.

I've made a simple implementation (it's not yet preserving pixel sizes, but does preserve tensor information), this is what it can do:

>>> import diplib as dip
>>> import numpy as np
>>> a = dip.ImageRead('examples/DIP.tif')
>>> a
<Color image (3x1 column vector, 3 elements, sRGB), UINT8, sizes {256, 256}>
>>> np.sum(a)
<Scalar image, UINT64, 0D>
>>> np.sum(a, axis=1)
<Color image (3x1 column vector, 3 elements, sRGB), UINT64, sizes {256}>
>>> np.sum(a, axis=2)
<Scalar image, UINT64, sizes {256, 256}>

Maybe I can get it to return the input object if it's a scalar, I don't yet know how to do this.

Implementation:

diplib/pydip/src/image.cpp

Lines 296 to 312 in 92071b2

	img.def( "__array_wrap__", []( dip::Image const& self, py::buffer& buf ) {
	dip::Image out = BufferToImage( buf, false );
	// Step 1: maybe one dimension is the tensor dimension
	dip::uint nTensor = self.TensorElements();
	for( dip::uint ii = 0; ii < out.Dimensionality(); ++ii ) {
	if( out.Size( ii ) == nTensor ) {
	out.SpatialToTensor( ii );
	out.ReshapeTensor( self.Tensor() );
	out.SetColorSpace( self.ColorSpace() );
	break;
	}
	}
	// Step 2: match dimensions by size (some might be missing, some might be new?)
	// TODO: how do we do this? It's an impossible task!
	// Step 3: copy over pixel sizes
	return out; // TODO: return a Python object so that we can return a scalar as not an image.
	});

@wcaarls What do you think of this? What behavior should it have? Maybe this is not even desirable?

[wcaarls]

Interesting! As an implementer I would not have immediately expected this behavior, but I could easily see a user wanting this. I'm a bit worried about the guessing of array dimensions, though. Would it be possible to use __array_ufunc__ to see which ufunc method was called, such that you know which dimensions were removed? The actual processing could be delegated to the original ufunc.

[crisluengo]

Yeah, I have been thinking about that. I'm guessing you'd need to know about all the ufuncs, and possibly check their parameters, to know what the output sizes will be. I don't think that is something we want to maintain.

I agree I don't like guessing at dimensions. A 3-channel image that has one dimension of size 3 might produce totally unexpected and wrong results with the code I wrote.

I'm not sure this can ever be done properly, but it was fun playing with for a while. :)

I actually don't mind the symmetry of the current situation, where calling a NumPy function on an image returns an array, and calling a DIPlib function on an array returns an image. Your "I would not have immediately expected this behavior" cements that. Thanks!

[crisluengo]

On the other hand, this will work:

a = dip.ImageRead('examples/trui')
b = dip.Image((256,), 1, "SFLOAT")
np.sum(a, out=np.asarray(b), axis=1)

BTW: I'm working to add bindings for the version of each function that has out as an argument rather than a return value. That works great too:

b = np.zeros((256,256), dtype=np.float64)
dip.Gauss(a, out=b, sigmas=1)

(And of course with a dip.Image object as output it works too.) It makes the out argument, and all following arguments, keyword arguments. By forcing writing out=..., it disambiguates the two function signatures.

It's a bit of work, and it adds a lot of duplication to the PyDIP source files, but it is very useful being able to specify an output buffer, or being able to some functions in-place.

[wcaarls]

Yeah, I have been thinking about that. I'm guessing you'd need to know about all the ufuncs, and possibly check their parameters, to know what the output sizes will be. I don't think that is something we want to maintain.

I'm not sure we need to know this on a per-ufunc basis. From a cursory glance, it seems their arguments that affect the output size are quite well-defined, particularly their methods.

[crisluengo]

Unrelated to the above. I have just found these two projects:

• Documentation generation for the Python bindings: https://github.com/pybind/pybind11_mkdoc

• More efficient alternative to pybind11 (faster compilation, faster execution, requires C++17, might not do everything we need it to do): https://github.com/wjakob/nanobind