Pillow (PIL) resize uses the NEAREST interpolation method by default.
Actually, the behavior differs between Pillow-3.3.3(October 5, 2016) and 3.4.0(October 4, 2016)

https://github.com/python-pillow/Pillow/releases?after=3.4.1

from PIL import Image
im = Image.open(sys.argv[1])
im = im.resize([2, 2], Image.NEAREST)
im.save("nearest.png")

A comparison of the results of resizing a 4x4 image to 2x2.

	image	dot enlargement
before resize
=< 3.3.3
3.4.0<=

why

Pillow's resize calls the function of affine transformation for scaling process.
Affine transformation is a convenient method to scale/rotate/translate, etc... all at once.

https://en.wikipedia.org/wiki/Affine_transformation

Scale	Rotate	Translate

As improvement of this affine transformation routine, when changing the calculation of pixel coordinates from the upper left corner of the pixel to calculation based on the pixel center, so it seems that the resize processing has been entangled. (I think it's too rough ^^;)

how

In coordinate calculation of which pixel of the original image to use, 0 (left or upper corner of pixel) start at Pillow-3.3.3 and earlier, 0.5(center of pixel) start at 3.4.0 and later.

<=3.3.3	3.4.0<=

When 0.5 is the pixel center, dividing one grid 0-0.999... into two divides it into 0-4.999... and 5.0-9.999..., in which case 0.5 is closer to the lower right. The difference in 3.4.0.

If it is at least 0.5-ε, what has not changed dramatically so far. (The problem comes out how to set ε value)

Pillow implementation

The simple code structure makes it easy to follow the process.

f_image => ImagingTransform with Affin
Affin is specified => ImagingTransformAffine
No rotation => ImagingScaleAffine

% tar tvfz Pillow-3.3.3.tar.gz
% cd Pillow-3.3.3

Starting from _imaging.c, concrete processing is inside libImaging.
_imaging.c: _resize

% grep ^_resize _imaging.c -A 32
_resize(ImagingObject* self, PyObject* args)
{
    Imaging imIn;
    Imaging imOut;

    int xsize, ysize;
    int filter = IMAGING_TRANSFORM_NEAREST;
    if (!PyArg_ParseTuple(args, "(ii)|i", &xsize, &ysize, &filter))
        return NULL;

    imIn = self->image;

    if (xsize < 1 || ysize < 1) {
        return ImagingError_ValueError("height and width must be > 0");
    }

    if (imIn->xsize == xsize && imIn->ysize == ysize) {
        imOut = ImagingCopy(imIn);
    }
    else if (filter == IMAGING_TRANSFORM_NEAREST) {
        double a[6];

        memset(a, 0, sizeof a);
        a[0] = (double) imIn->xsize / xsize;
        a[4] = (double) imIn->ysize / ysize;

        imOut = ImagingNew(imIn->mode, xsize, ysize);

       imOut = ImagingTransform(
            imOut, imIn, IMAGING_TRANSFORM_AFFINE,
            0, 0, xsize, ysize,
            a, filter, 1);

libImaging/Geometry.c: ImagingTransform


% grep "^ImagingTransform(" libImaging/Geometry.c  -A 10
ImagingTransform(Imaging imOut, Imaging imIn, int method,
                 int x0, int y0, int x1, int y1,
                 double a[8], int filterid, int fill)
{
    ImagingTransformMap transform;

    switch(method) {
    case IMAGING_TRANSFORM_AFFINE:
        return ImagingTransformAffine(
            imOut, imIn, x0, y0, x1, y1, a, filterid, fill);
        break;

libImaging/Geometry.c: ImagingTransformAffine

% grep ^ImagingTransformAffine libImaging/Geometry.c  -A 25
ImagingTransformAffine(Imaging imOut, Imaging imIn,
                       int x0, int y0, int x1, int y1,
                       double a[6], int filterid, int fill)
{
    /* affine transform, nearest neighbour resampling, floating point
       arithmetics*/

    ImagingSectionCookie cookie;
    int x, y;
    int xin, yin;
    int xsize, ysize;
    double xx, yy;
    double xo, yo;

    if (filterid || imIn->type == IMAGING_TYPE_SPECIAL) {
        return ImagingGenericTransform(
            imOut, imIn,
            x0, y0, x1, y1,
            affine_transform, a,
            filterid, fill);
    }

    if (a[1] == 0 && a[3] == 0)
        /* Scaling */
        return ImagingScaleAffine(imOut, imIn, x0, y0, x1, y1, a, fill);

libImaging/Geometry.c:ImagingScaleAffine


% grep ^ImagingScaleAffine libImaging/Geometry.c  -A 50
ImagingScaleAffine(Imaging imOut, Imaging imIn,
                   int x0, int y0, int x1, int y1,
                   double a[6], int fill)
{
    /* scale, nearest neighbour resampling */

    ImagingSectionCookie cookie;
    int x, y;
    int xin;
    double xo, yo;
    int xmin, xmax;
    int *xintab;

    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
        return (Imaging) ImagingError_ModeError();

    ImagingCopyInfo(imOut, imIn);

    if (x0 < 0)
        x0 = 0;
    if (y0 < 0)
        y0 = 0;
    if (x1 > imOut->xsize)
        x1 = imOut->xsize;
    if (y1 > imOut->ysize)
        y1 = imOut->ysize;

    /* malloc check ok, uses calloc for overflow */
    xintab = (int*) calloc(imOut->xsize, sizeof(int));
    if (!xintab) {
        ImagingDelete(imOut);
        return (Imaging) ImagingError_MemoryError();
    }

    xo = a[2];
    yo = a[5];

    xmin = x1;
    xmax = x0;

    /* Pretabulate horizontal pixel positions */
    for (x = x0; x < x1; x++) {
        xin = COORD(xo);
        if (xin >= 0 && xin < (int) imIn->xsize) {
            xmax = x+1;
            if (x < xmin)
                xmin = x;
            xintab[x] = xin;
        }
        xo += a[0];
    }

The initial value of this xo and yo changes with Pillow-3.4.0.

% tar tvfz Pillow-3.3.4.tar.gz
% cd Pillow-3.3.4

% grep ^ImagingScaleAffine libImaging/Geometry.c  -A 50 | grep "xo =" -A 4
    xo = a[2] + a[0] * 0.5;
    yo = a[5] + a[4] * 0.5;

    xmin = x1;
    xmax = x0;

finaly

resize We explained the situation about the default method NEAREST, but BILINEAR, BICUBIC, LANCZOS etc. also have differences due to the same reason. be careful.

reference