First, a big thank you for taking this on. This is an important part of matplotlib that needs improvement.

One concern: This is effectively decimating (i.e. ignoring data), not downsampling (i.e. averaging across neighboring values), correct? I understand that might be acceptable in some domains, but probably not in others (particularly in astronomy I know this would almost never be the right thing). That's not a reason to exclude it, but it should be documented clearly that that is what this feature does.

The bigger problem matplotlib has currently is that images are processed by:

• converting floating point values to RGB (colormapping)
• interpolating in RGB space

It's done this way historically because we're using Agg, which has a bunch of interpolation algorithms in place that were easy to grab and use, and these only supported RGB. With the latest SVN version of Agg, this is reportedly no longer the case, but I spent some time a few months ago looking at this and it seems to not work as advertised.

But assuming someone can do the work so instead we do:

• interpolate in floating point
• colormap floating point -> RGB

then what you've done here just becomes one of many interpolation algorithms (and probably the fastest one).

Are you interested in helping with this?

See latest commit, which allows the choice of downsampling method. The commit message includes an example where decimation vs. e.g. averaging produces significantly different results. (I wrote the decimating code first mostly for simplicity (with speed as a secondary benefit) but index_tricks comes to the rescue as usual.)
I think this issue is pretty different from interpolation (and I certainly don't intend to rewrite all of agg's interpolation methods :-)), here instead we're grouping multiple values into a single one.

I think long term we really want to do proper interpolation on the input (greyscale) image, and this is just a half measure. I've already solved the C++ side of it in matplotlib/aggravate, now we just need to make the Python side use it. I think that will ultimately be more flexible than this and also be the same speed when interpolation == nearest. Would you like to help with that effort?

To be honest I don't care much, or at least have very little intuition, about interpolation (I only believe in interpolation="none" :-)) but more importantly I still fail to see the connection between the two cases -- other than in the case where you have too many data points per pixel in one direction (you need to downsample there) and not enough in the other (you need to interpolate there) (but then you can always do the two parts sequentially).
A quick look (... limited by the quality of the docs but heh) at Agg's resampling algorithms suggests that they do not support downsampling, only upsampling, right?

A quick look (... limited by the quality of the docs but heh) at Agg's resampling algorithms suggests that they do not support downsampling, only upsampling, right?

Nope, they address both. It's really just generic "resampling" where each pixel in the output is projected onto the input, and the average of that rectangle taken (modulo the filtering kernel).

But do we always want to use averaging? I'd say this is pretty application dependent and certainly at least averaging, geometric averaging, median and possibly decimation all have legitimate uses there.

How does this connect with the new image infrastructure in #5718

@mdboom suggests that #5718 handles this issue as a special case, but I haven't really followed what is happening on that issue.

This is a special case of #5718 (see my pathological example).

This is basically a performance improvement for plotting very large images through decimation (which only really does the right thing if interpolation=='nearest').

To test, I'm using the following:

from matplotlib import pyplot as plt
import numpy as np

Z = np.random.rand(10000, 10000)
plt.imshow(Z, interpolation='nearest')
plt.savefig("test.png")

(If you want to reproduce the numbers, make sure you run it twice and ignore the first run to remove the effects of regenerating the font cache).

I get:

master: 7.24s
#5718: 3.61s
#5602: 2.22s

So, #5718 is much better than master already, by virtue of downsampling first, and then colormapping something much smaller. I'm a bit surprised by the remaining advantage of this over #5718, since it should boil down to effectively the same thing computationally... I'll need to do some further investigation. I'd still like to find a more general solution than this if possible.

As best I can tell from profiling, the difference in time is mostly in norm. This PR runs norm on the decimated data, which unfortunately isn't correct, since it could miss local maxima and minima in the data that way if they fall in the cracks. I think we have to norm (or at least find min and max) over the entire original data, which obviously takes considerably more memory bandwidth which is I believe what accounts for the difference in timing.

Two quick notes:

• In its current state, the patch allows "fairly easily" other downsampling methods than decimation:

                # This is where downsampling happens.
                A = A[:, :, 0, 0] # or A = A.mean(axis=(-1, -2)), or etc.

(image.py line 239-240). Not that it matters because of the following point:

• You are correct that normalization needs to happen on the whole image before resampling (typically: if you use a Power norm, the proper way to average your values is probably a geometric mean, so you may as well normalize to log scale and do an arithmetic mean after). So only the decimation method can benefit from normalization-after-resampling. (Also, I guess the difference in timing isn't just from computing the min and max, but also from doing the normalization itself on the whole vs the downsampled array.)

@mdboom, @anntzer, I'm not following this in detail, but here are a couple thoughts that might be relevant to the discussion:

1. The computation of min and max is going through ma.min and ma.max, which probably has two sources of inefficiency: first, ma operations tend to be slow for a variety of reasons; and second, because it is unnecessarily making two passes instead of calculating min and max in a single pass. It's unfortunate that numpy has ptp(), but not a function that returns min and max from a single pass. We could use our own C++ code for this (handling the possible mask at the same time).
2. I haven't checked, but in almost anything related to color handling in floating point, if we are using double precision we are wasting memory. I think there is a place where we push things toward single precision, but I don't know whether we are enforcing single precision.
   Sorry not to be able to look into this more closely for now.
   Eric

1. The computation of min and max is going through ma.min and ma.max, which probably has two sources of inefficiency: first, ma operations tend to be slow for a variety of reasons; and second, because it is unnecessarily making two passes instead of calculating min and max in a single pass. It's unfortunate that numpy has ptp(), but not a function that returns min and max from a single pass. We could use our own C++ code for this (handling the possible mask at the same time).

Yeah -- this seems like something that Numpy really should grow. I wonder if they'd be amenable to a pull request on that. We'd still have to have our own answer in the meantime, though.

2. I haven't checked, but in almost anything related to color handling in floating point, if we are using double precision we are wasting memory. I think there is a place where we push things toward single precision, but I don't know whether we are enforcing single precision.

In the common case of single plane data that is color mapped, we color map to 8-bit ints (the colormap code can also output to floats, but in practice that code path is rarely used). #5718 improves things further when the input data is already color: in that case, if the input data is 8-bit int, we no longer require passing through float, but can preserve it as 8-bit all the way down.

There have been talk of a np.minmax() function for years, but it has never
moved beyond that. I think they would be amendable to such a PR.

Something else that is even worse about np.ma.min() and such... they make
copies of the array (sans masked parts)! I was certain there was a github
issue filed for that, but I can't find it at the moment.

On Tue, Jan 5, 2016 at 2:11 PM, Michael Droettboom <notifications@github.com

wrote:

1. The computation of min and max is going through ma.min and ma.max,
   which probably has two sources of inefficiency: first, ma operations tend
   to be slow for a variety of reasons; and second, because it is
   unnecessarily making two passes instead of calculating min and max in a
   single pass. It's unfortunate that numpy has ptp(), but not a function that
   returns min and max from a single pass. We could use our own C++ code for
   this (handling the possible mask at the same time).

Yeah -- this seems like something that Numpy really should grow. I wonder
if they'd be amenable to a pull request on that. We'd still have to have
our own answer in the meantime, though.

2. I haven't checked, but in almost anything related to color handling in
   floating point, if we are using double precision we are wasting memory. I
   think there is a place where we push things toward single precision, but I
   don't know whether we are enforcing single precision.

In the common case of single plane data that is color mapped, we color map
to 8-bit ints (the colormap code can also output to floats, but in practice
that code path is rarely used). #5718
#5718 improves things
further when the input data is already color: in that case, if the input
data is 8-bit int, we no longer require passing through float, but can
preserve it as 8-bit all the way down.

—
Reply to this email directly or view it on GitHub
#5602 (comment)
.

The following shows where some major slowdowns are coming from the masked array support. It certainly would be nice to get all that sped up within numpy...

In [1]: import numpy as np
In [2]: xx = np.random.randn(10000, 10000)
In [3]: timeit(xx.min())
10 loops, best of 3: 47 ms per loop
In [4]: timeit(xx.ptp())
10 loops, best of 3: 94.5 ms per loop
In [7]: yy = np.ma.masked_greater(xx, 2)
In [8]: timeit(yy.min())
1 loops, best of 3: 508 ms per loop
In [9]: timeit np.nanmin(yy.filled(np.nan))
1 loops, best of 3: 577 ms per loop
In [10]: timeit yy.filled(np.nan)
1 loops, best of 3: 457 ms per loop
In [12]: timeit np.ma.max(xx)
10 loops, best of 3: 47.1 ms per loop

It looks like the implementation of ptp() is not very efficient.
Masked arrays are slow only when there are masked values, so it doesn't look like min and max are responsible for the slow norming in your test case--correct?

Masked arrays are slow only when there are masked values, so it doesn't look like min and max are responsible for the slow norming in your test case--correct?

Yeah -- in my test case, there's nothing masked out...

(Also, I guess the difference in timing isn't just from computing the min and max, but also from doing the normalization itself on the whole vs the downsampled array.)

Yes -- there may be a way to take advantage of that in the context of #5718. The tricky bit is that normalization happens in Python/Numpy (and through a nice plug-in API where different normalization functions can be written and interchanged), and resampling happens in C++. Plugging normalization into the resampling step might be tricky, though. #5718 already cuts the run time in half -- we know from this PR that there's more that could probably be done wrt to speed, but I'm not certain we should let that hold up #5718.

How does the new image machinery play with this? I have noticed it seems to perform quite well even with very large images. Should this PR be closed?

Thanks for pinging this. Indeed, the image refactor was in part a response to this -- to try to do this faster without resorting to decimation. I think we should probably close this unless further testing reveals new problems with the new code.

Reopening per the discussion at https://gitter.im/matplotlib/matplotlib?at=5c92437e49fdaa0d81e71ec9.
(Obviously the PR needs a big update, this is just to track the issue.)

Image taken from http://www.imagemagick.org/Usage/filter/nicolas/rings_sm_orig.gif The article is also quite useful: http://www.imagemagick.org/Usage/filter/nicolas/

import matplotlib.pyplot as plt
import numpy as np

a = plt.imread('rings_sm_orig.gif')

print(np.shape(a))
# 200x200
fig, ax = plt.subplots(figsize=(4, 4))
ax.imshow(a, cmap='gray')
ax.set_position([0, 0, 1, 1])
fig.savefig('dpi50.png', dpi=50)
fig.savefig('dpi30.png', dpi=30)


fig, ax = plt.subplots(figsize=(4, 4))
ax.imshow(a, cmap='gray', interpolation='lanczos', filterrad=(2*50/30))
ax.set_position([0, 0, 1, 1])
fig.savefig('filtdpi50.png', dpi=50)
fig.savefig('filtdpi30.png', dpi=30)

dpi50.png:

dpi30.png

0.png

With filter:

filtdpi50.png

filtdpi30.png

Conclusion:

Agg seems to do the right thing if we give it the appropriate filter. The question is if we apply this filter automatically if we know we are downsampling....

The issue is better described in #13724 now, so I'll close this.