Incorrect upsampling with Image.BOX

[99991]

What did you do?

I created a white image and resized it using Image.BOX resampling.

What did you expect to happen?

I expected a white image:

What actually happened?

But instead, I got a white image with vertical black lines:

What are your OS, Python and Pillow versions?

• OS: Debian 10
• Python: 3.7.3
• Pillow: 6.2.1

from PIL import Image

# create a white image
image = Image.new('RGB', (8, 100), (255, 255, 255))
# resize image
image = image.resize((100, 100), Image.BOX)
# show white image with black vertical lines
image.show()

I think a good starting point to find this bug would be to check if any of the horizontal filtering coefficients are zero.

[homm]

I can confirm. The problem is that BOX, unlike other filters, is a discontinuous function. It equals 1 in range (-0.5, 0.5) and equals 0 on the rest range. What is it equal to in -0.5 and 0.5? I don't know exactly, because it should be exactly 1 long, so if it equals 1, it will be a bit longer than 1, and if it will be 0, it will be slightly shorter than 1.

Current implementation defines box filter as 1 on [-0.5, 0.5) range (0.5 is not included).

static inline double box_filter(double x)
{
    if (x >= -0.5 && x < 0.5)
        return 1.0;
    return 0.0;
}

I'll try to find the way how to fix this gently.

[99991]

For reference, OpenCV's INTER_AREA filter builds the weighting table like this.

[homm]

Oh my god, I believe I found a fix and can mathematically prove it.

The fix is very simple. Instead of defining the BOX filter as 1 on [-0.5, 0.5), define it 1 on (-0.5, 0.5].

static inline double box_filter(double x)
{
    if (x > -0.5 && x <= 0.5)
        return 1.0;
    return 0.0;
}

The main question is "how we can be sure that 0.0 is the correct value for box_filter(-0.5) while it wasn't for box_filter(0.5)?".

To prove this I'm going to prove that x argument for box_filter could never be -0.5.

I will use the following lines from the code:

    filterscale = max(1.0, scale)
    support = filterp->support * filterscale;
    ss = 1.0 / filterscale;
    xmin = (int) (center - support + 0.5);
    for (x = 0; x < xmax; x++) {
        double w = filterp->filter((x + xmin - center + 0.5) * ss);
        k[x] = w;
        ww += w;
    }

Start

For x == -0.5 in box_filter function, this equation should be true:
(x + xmin - center + 0.5) * ss == -0.5

ss is 1.0 / filterscale, where filterscale is max(1.0, scale):
x + xmin - center + 0.5 == -0.5 * max(1.0, scale)

Now we have two branches, first with 0.0 < scale <= 1.0 and filterscale == 1.0, second with
scale > 1.0 and filterscale == scale

0.0 < scale <= 1.0 case

x + xmin - center + 0.5 == -0.5

xmin is defined as int(center - support + 0.5):
x + int(center - support + 0.5) - center + 1 == 0

support is filterp->support (which is 0.5 for BOX filter) multiplied by filterscale, which is 1.0:
x + int(center - 0.5 + 0.5) - center + 1 == 0

Rearrange:
x + 1 == center - int(center)

smth - int(smth) is a fraction part of smth. In the end:
x + 1 == frac(center)

scale > 1.0 case

x + xmin - center + 0.5 == -0.5 * scale

xmin is defined as int(center - support + 0.5):
x + int(center - support + 0.5) - center + 0.5 == -0.5 * scale

support is filterp->support (which is 0.5 for BOX filter) multiplied by filterscale, which is scale:
x + int(center - 0.5 * scale + 0.5) - center + 0.5 == -0.5 * scale

Rearrange:
x + 1 == center - 0.5 * scale + 0.5 - int(center - 0.5 * scale + 0.5)

So, we came up to the same equation:
x + 1 == frac(center - 0.5 * scale + 0.5)

Summary

We have two cases:
x + 1 == frac(center)
x + 1 == frac(center - 0.5 * scale + 0.5)

Since x is some positive integer, and a fraction is always in [0.0, 1.0) range, but newer equals 1.0, both equations are impossible.

Therefore x could newer be -0.5 in box_filter function.