Axes.imshow draws invalid color at value is 0 when max of 'X' not equal to vmax

Bug report

Bug summary

Axes.imshow draws invalid color at cells which value is 0 when max value of parameter 'X' is not equal to parameter 'vmax'

Code for reproduction

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import numpy as np

cmap = plt.cm.get_cmap('autumn_r')
cmap.set_under(color='lightgrey')

fig = plt.figure()
ax = fig.add_subplot()
ax.grid(which='major', color='white', linestyle='dotted', alpha=0.7)
ax.grid(which='minor', color='white', linestyle='dotted', alpha=0.5)

data = np.array([[-1, -1, -1, 0, 0, 0, 0, 43, 79, 95, 66, 1, -1, -1, -1, 0, 0, 0, 34]])
# data = np.array([[-1, -1, -1, 0, 0, 0, 0, 43, 79, 100, 66, 1, -1, -1, -1, 0, 0, 0, 34]])
im = ax.imshow(data, aspect='auto', cmap=cmap, vmin=0, vmax=100)
ax.xaxis.set_major_locator(ticker.MultipleLocator())
ax.set_yticks([0])
fig.colorbar(im)

plt.show()

Actual outcome

I expect yellow at cells 3, 4, 5, 6, 15, 16, 17. but they are grey.

Expected outcome

if using commented data(max(data) == vmax), got output below. that is valid.

Matplotlib version

Operating system: Windows10
Matplotlib version: 3.1.3
Matplotlib backend (print(matplotlib.get_backend())): TkAgg
Python version: 3.7.4

I assumed you meant min value not equal to vmin.

Here you explicitly requested that values below vmin are in lightgrey (`cmap.set_under('lightgrey')), so I think this is behaving as expected?

I think the report is very consistent. But we can make it more clear: Compare how the same data is plotted via imshow and pcolormesh. Of course we would expect both images to be equal.

import matplotlib.pyplot as plt
import numpy as np

cmap = plt.cm.get_cmap('autumn_r')
cmap.set_under(color='lightgrey')

norm = plt.Normalize(0, 100)

data = np.array([[-1, -1, -1, 0, 0, 0, 95, 95, 95]])

fig, (ax1, ax2) = plt.subplots(ncols=2)

im1 = ax1.imshow(data, aspect='auto', cmap=cmap, norm=norm)
im2 = ax2.pcolormesh(data, cmap=cmap, norm=norm)

fig.colorbar(im1, ax=ax1)
fig.colorbar(im2, ax=ax2)

plt.show()

This works with matplotlib 2.0.2, but fails with matplotlib 2.1.2 or higher. (Unfortunately, I'm unable to bisect)

Will mark as release critical for 3.3, though one might even consider this a bad enough bug to fix earlier.

My knee-jerk reaction is that this a numerical stability issue.

Possibly related to #10072 / #10133 . Chasing through the linked bug reports and PRs is a good place to start. b503d6a (2.1.x branch) or b511bb2 (on the master branch) are likely the offending commits.

That this now "works" for pcolormesh and not imshow suggests that the other bug reported in #10072 is still present in pcolormesh.

Bisect says 12c27f3, i.e., #8966.

@anntzer

I think the problem is

matplotlib/lib/matplotlib/image.py

Lines 469 to 484 in 15c0d7c

    
           if a_min != a_max: 
        
               A_scaled /= ((a_max - a_min) / 0.8) 
        
           A_scaled += 0.1 
        
           # resample the input data to the correct resolution and shape 
        
           A_resampled = _resample(self, A_scaled, out_shape, t) 
        
           # done with A_scaled now, remove from namespace to be sure! 
        
           del A_scaled 
        
           # un-scale the resampled data to approximately the 
        
           # original range things that interpolated to above / 
        
           # below the original min/max will still be above / 
        
           # below, but possibly clipped in the case of higher order 
        
           # interpolation + drastically changing data. 
        
           A_resampled -= 0.1 
        
           if a_min != a_max: 
        
               A_resampled *= ((a_max - a_min) / 0.8) 
        
           A_resampled += a_min

The issue that code is trying to solve is clipping at the edges of the interpolation range so we take the input data, scale it to [.1, .9] , run it through the re-sampling routine, re sample (up or down) the data to match the screen resolution, then shift it back to the input range, then run it through the normalization (linear, log, etc), then color map it from [0, 1] -> RGBA, then draw that to the screen.

Pulling out the scale down / scale up round tripping

import numpy as np
import matplotlib.pyplot as plt

def test_broken(top, scaled_dtype=np.float32):
    A = np.array([-1, 0, top], dtype=float)
    a_min = A.min()
    a_max = A.max()
    A_scaled = np.array(A, dtype=scaled_dtype)
    A_scaled -= a_min

    A_scaled /= (a_max - a_min) / 0.8
    A_scaled += 0.1

    # un-scale the resampled data to approximately the
    # original range things that interpolated to above /
    # below the original min/max will still be above /
    # below, but possibly clipped in the case of higher order
    # interpolation + drastically changing data.
    A_scaled = np.array(A_scaled)
    A_scaled -= 0.1
    A_scaled *= (a_max - a_min) / 0.8
    A_scaled += a_min

    return (A_scaled - A)[1]




fig, (ax1, ax2) = plt.subplots(2, constrained_layout=True)

for ax, sdt in zip((ax1, ax2), (np.float32, np.float64)):
    data = [test_broken(n, scaled_dtype=sdt) for n in range(250)]
    ax.plot(data, label=str(sdt))
    ax.set_title(f'round-trip error of 0 with {sdt!r}')
    ax.set_ylabel('error')
    ax.set_xlabel('"0" value')

plt.show()

we see that we systematically miss the round trip of the 0 in the array (both above and below) based on what the maximum limits are (it would be interesting to extend this to do the matrix of varying the low value as well).

This looks like it is fundamental issue with float precision so tweaking the number on the way we do the scaling (the .1 and .8 constants could be tweaked) will avoid the problem in all cases. I am not enough of an expert on numerical computation to know if there are things we could be doing in this code to reduce the error.

We can not do the re-scaling for interpolation clipped to the vmin/vmax range because that causes very out-of-range points to contribute less (the vmid logic in the code is a compromise to handle very small vmin/vmax on a very large background).

I am going to move this to 3.4 as I suspect the real solution here is to (as @anntzer has been advocating for a while) is replace the Agg based interpolation with something that does not require us to pre-scale to [0, 1].

The very blunt work-around, if you know your data is going to be integers is to set the vmin/vmax to be a float just below and just above the min / max values.

🤦 Sorry, lost track of what I was doing, your right, it does not eliminate the problem.

So maybe something like this, with a smaller buffer, and then in colors.py around line 571:

        eps = 1e-16  * self.N   
        xa[xa > self.N - 1 + eps] = self._i_over
        xa[xa < -eps] = self._i_under

though there is a funky statement above on line 561 that also needs to be changed.

EDIT: where I'm not sure 1e-16 is the right number at all, and I guess from the above its not for 32-bit images.

Maybe the right thing to do is to run vmin and vmax through the same pipeline?

I don't know if that helps? I guess that helps for the exact x=vmin case, but not sure if it helps for x = vmin + eps case, which I think could theoretically end up on the wrong side of vmin?

I hope that if a < b is true is the input data, then the worst case scenario after the transform is a <= b (and a > b should be impossible). Hence if x = vmin + eps it should still stay on "the right side" of the mapping?

Good question, but afraid I don't know the answer!

Actually, I'm vaguely confused now. Why do we have to rescale to 0.1 0.9 rather than e.g. 0/1? We are separately resampling the mask anyways, so we should be able to use that info anyways to know what to mask out at the end, no? Also, what's the way resampling between an non-finite and a finite value is handled?

Lets just say vmin=-1 and vmax=1. If we renormalize from 0-1, then a value at -1.1 gets set to 0 and we lose the idea that it is "under" -1. But we want to keep track of the fact that its under, so hence the kludge.

MNT: also transform vmin/vmax when rendering imshow

But what where we doing with nans? They were mapped to something, no? And then masked out later, because we also resample the nan mask. So why can't we have an under/over mask as well?

Experimentally, running the same rescaling logic on vmin/vmax and the temporarily adjusting the vmin/vmax on the norm does indeed fix the OP.

We did have an over/under mask at one point, but we had to back out of that because when you use the non-trivial interpolations you may get values that are with in the allowed gamut but include values from the original data with is out of gamut. At one point we did mask any pixel which was interpolated from an out-of-gamut input as out of range and it very quickly got bug reports so we have to do the over/under classification after we do the resampling. We do poison invalid (nan or inf) points.

For completeness

import numpy as np
import matplotlib.pyplot as plt

exp = 5
offset = 1/(2**exp)
frac = 1 - 2 * offset


def test_broken(top, frac, offset, scaled_dtype=np.float32):
    A = np.array([-1, 0, top], dtype=float)
    a_min = A.min()
    a_max = A.max()
    A_scaled = np.array(A, dtype=scaled_dtype)
    A_scaled -= a_min

    A_scaled /= (a_max - a_min) / frac
    A_scaled += offset

    # un-scale the resampled data to approximately the
    # original range things that interpolated to above /
    # below the original min/max will still be above /
    # below, but possibly clipped in the case of higher order
    # interpolation + drastically changing data.
    A_scaled = np.array(A_scaled)
    A_scaled -= offset
    A_scaled *= (a_max - a_min) / frac
    A_scaled += a_min

    return (A_scaled - A)[1]




fig, (ax1, ax2) = plt.subplots(2, constrained_layout=True)

for ax, sdt in zip((ax1, ax2), (np.float32, np.float64)):
    for exp in [2, 3, 4, 5]:
        offset = 1/(2**exp)
        frac = 1 - 2 * offset
        data = [test_broken(n, frac, offset, scaled_dtype=sdt) for n in range(1000)]
        ax.plot(data, label=f'{offset}')
        ax.set_title(f'round-trip error of 0 with {sdt!r}')
        ax.set_ylabel('error')
        ax.set_xlabel('top value')
    ax.legend(ncol=2)
plt.show()

is the corrected code for the last set of images I posted.

Navigation Menu

Repository navigation

Axes.imshow draws invalid color at value is 0 when max of 'X' not equal to vmax #16910

Bug report

12 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	if a_min != a_max:
	A_scaled /= ((a_max - a_min) / 0.8)
	A_scaled += 0.1
	# resample the input data to the correct resolution and shape
	A_resampled = _resample(self, A_scaled, out_shape, t)
	# done with A_scaled now, remove from namespace to be sure!
	del A_scaled
	# un-scale the resampled data to approximately the
	# original range things that interpolated to above /
	# below the original min/max will still be above /
	# below, but possibly clipped in the case of higher order
	# interpolation + drastically changing data.
	A_resampled -= 0.1
	if a_min != a_max:
	A_resampled *= ((a_max - a_min) / 0.8)
	A_resampled += a_min

Axes.imshow draws invalid color at value is 0 when max of 'X' not equal to vmax #16910

Description

Bug report

Activity

anntzer commented on Mar 26, 2020

ImportanceOfBeingErnest commented on Mar 26, 2020

ImportanceOfBeingErnest commented on Mar 26, 2020

tacaswell commented on Mar 26, 2020

tacaswell commented on Mar 26, 2020

QuLogic commented on Apr 29, 2020

tacaswell commented on Jun 12, 2020

12 remaining items

tacaswell commented on Jun 15, 2020

jklymak commented on Jun 15, 2020

tacaswell commented on Jun 15, 2020

jklymak commented on Jun 15, 2020

tacaswell commented on Jun 15, 2020

jklymak commented on Jun 15, 2020

anntzer commented on Jun 15, 2020

jklymak commented on Jun 15, 2020

anntzer commented on Jun 15, 2020

tacaswell commented on Jun 15, 2020

tacaswell commented on Oct 7, 2020

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions