Pillow(PIL) の resize はデフォルトで NEAREST 補間方式を使います。
これ、実は Pillow-3.3.3 (2016年10月5日)以前と、3.4.0(2016年10月4日)以降で動きが異なります。

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")

4x4 画像を 2x2 にリサイズした結果の比較です。

	画像	(ドット拡大)	対応するピクセル
リサイズ前画像
3.3.3以前
3.4.0以降

why (何故？)

Pillow のリサイズはアフィン変換の関数を呼んで拡大縮小処理を任せます。
アフィン変換は拡大縮小、回転、移動などを一度に行える便利な手法です。

https://ja.wikipedia.org/wiki/アフィン写像

Scale	Rotate	Translate

このアフィン変換の改善としてピクセル座標の計算がピクセル左上隅だったのをピクセル中心ベースに計算するよう変更した際、resize 処理はその巻き添えになったようです。(正直、雑すぎると思います ^^;)

how (どのように？)

元画像のどのピクセルを使うかの座標計算で、Pillow-3.3.3 までは 0 (ピクセルの左or上隅)スタートだったのが、3.4.0 以降では 0.5 開始にしています。

3.3.3 以前	3.4.0 以降

0.5 をピクセル中心とする場合、一つのグリッド 0〜0.999... を2分割すると 0〜4.999... と 5.0〜9.999... に分かれ、どちらかというと 0.5 は右下寄りになります。その違いが 3.4.0 での差異になります。
せめて 0.5 - ε にしてくれれば、ここまで劇的に変わらなかったものを。(ε 値をいくつにするか問題が出ててくるけど)

Pillow 実装

素直な構造で、処理を追いやすいです。

_image から Affin 指定で ImagingTransform
Affin 指定なので ImagingTransformAffine
回転が無いので ImagingScaleAffine

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

_imaging.c が起点で、具体的な処理は 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];
    }

この xo と yo の初期値が Pillow-3.4.0 で変わります。

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

% 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;

最後に

resize デフォルトの NEAREST について状況を説明しましたが、BILINEAR, BICUBIC, LANCZOS 等も同様の理由により違いが生じます。ご注意下さい。

参考