diff --git a/Tests/images/radial_gradients.png b/Tests/images/radial_gradients.png
new file mode 100644
index 000000000..39a02fbbf
Binary files /dev/null and b/Tests/images/radial_gradients.png differ
diff --git a/Tests/test_image_reduce.py b/Tests/test_image_reduce.py
new file mode 100644
index 000000000..d8f2ce1ec
--- /dev/null
+++ b/Tests/test_image_reduce.py
@@ -0,0 +1,244 @@
+from PIL import Image, ImageMath, ImageMode
+
+from .helper import PillowTestCase, convert_to_comparable
+
+
+class TestImageReduce(PillowTestCase):
+    # There are several internal implementations
+    remarkable_factors = [
+        # special implementations
+        1,
+        2,
+        3,
+        4,
+        5,
+        6,
+        # 1xN implementation
+        (1, 2),
+        (1, 3),
+        (1, 4),
+        (1, 7),
+        # Nx1 implementation
+        (2, 1),
+        (3, 1),
+        (4, 1),
+        (7, 1),
+        # general implementation with different paths
+        (4, 6),
+        (5, 6),
+        (4, 7),
+        (5, 7),
+        (19, 17),
+    ]
+
+    @classmethod
+    def setUpClass(cls):
+        cls.gradients_image = Image.open("Tests/images/radial_gradients.png")
+        cls.gradients_image.load()
+
+    def test_args_factor(self):
+        im = Image.new("L", (10, 10))
+
+        self.assertEqual((4, 4), im.reduce(3).size)
+        self.assertEqual((4, 10), im.reduce((3, 1)).size)
+        self.assertEqual((10, 4), im.reduce((1, 3)).size)
+
+        with self.assertRaises(ValueError):
+            im.reduce(0)
+        with self.assertRaises(TypeError):
+            im.reduce(2.0)
+        with self.assertRaises(ValueError):
+            im.reduce((0, 10))
+
+    def test_args_box(self):
+        im = Image.new("L", (10, 10))
+
+        self.assertEqual((5, 5), im.reduce(2, (0, 0, 10, 10)).size)
+        self.assertEqual((1, 1), im.reduce(2, (5, 5, 6, 6)).size)
+
+        with self.assertRaises(TypeError):
+            im.reduce(2, "stri")
+        with self.assertRaises(TypeError):
+            im.reduce(2, 2)
+        with self.assertRaises(ValueError):
+            im.reduce(2, (0, 0, 11, 10))
+        with self.assertRaises(ValueError):
+            im.reduce(2, (0, 0, 10, 11))
+        with self.assertRaises(ValueError):
+            im.reduce(2, (-1, 0, 10, 10))
+        with self.assertRaises(ValueError):
+            im.reduce(2, (0, -1, 10, 10))
+        with self.assertRaises(ValueError):
+            im.reduce(2, (0, 5, 10, 5))
+        with self.assertRaises(ValueError):
+            im.reduce(2, (5, 0, 5, 10))
+
+    def test_unsupported_modes(self):
+        im = Image.new("P", (10, 10))
+        with self.assertRaises(ValueError):
+            im.reduce(3)
+
+        im = Image.new("1", (10, 10))
+        with self.assertRaises(ValueError):
+            im.reduce(3)
+
+        im = Image.new("I;16", (10, 10))
+        with self.assertRaises(ValueError):
+            im.reduce(3)
+
+    def get_image(self, mode):
+        mode_info = ImageMode.getmode(mode)
+        if mode_info.basetype == "L":
+            bands = [self.gradients_image]
+            for _ in mode_info.bands[1:]:
+                # rotate previous image
+                band = bands[-1].transpose(Image.ROTATE_90)
+                bands.append(band)
+            # Correct alpha channel by transforming completely transparent pixels.
+            # Low alpha values also emphasize error after alpha multiplication.
+            if mode.endswith("A"):
+                bands[-1] = bands[-1].point(lambda x: int(85 + x / 1.5))
+            im = Image.merge(mode, bands)
+        else:
+            assert len(mode_info.bands) == 1
+            im = self.gradients_image.convert(mode)
+        # change the height to make a not-square image
+        return im.crop((0, 0, im.width, im.height - 5))
+
+    def compare_reduce_with_box(self, im, factor):
+        box = (11, 13, 146, 164)
+        reduced = im.reduce(factor, box=box)
+        reference = im.crop(box).reduce(factor)
+        self.assertEqual(reduced, reference)
+
+    def compare_reduce_with_reference(self, im, factor, average_diff=0.4, max_diff=1):
+        """Image.reduce() should look very similar to Image.resize(BOX).
+
+        A reference image is compiled from a large source area
+        and possible last column and last row.
+        +-----------+
+        |..........c|
+        |..........c|
+        |..........c|
+        |rrrrrrrrrrp|
+        +-----------+
+        """
+        reduced = im.reduce(factor)
+
+        if not isinstance(factor, (list, tuple)):
+            factor = (factor, factor)
+
+        reference = Image.new(im.mode, reduced.size)
+        area_size = (im.size[0] // factor[0], im.size[1] // factor[1])
+        area_box = (0, 0, area_size[0] * factor[0], area_size[1] * factor[1])
+        area = im.resize(area_size, Image.BOX, area_box)
+        reference.paste(area, (0, 0))
+
+        if area_size[0] < reduced.size[0]:
+            self.assertEqual(reduced.size[0] - area_size[0], 1)
+            last_column_box = (area_box[2], 0, im.size[0], area_box[3])
+            last_column = im.resize((1, area_size[1]), Image.BOX, last_column_box)
+            reference.paste(last_column, (area_size[0], 0))
+
+        if area_size[1] < reduced.size[1]:
+            self.assertEqual(reduced.size[1] - area_size[1], 1)
+            last_row_box = (0, area_box[3], area_box[2], im.size[1])
+            last_row = im.resize((area_size[0], 1), Image.BOX, last_row_box)
+            reference.paste(last_row, (0, area_size[1]))
+
+        if area_size[0] < reduced.size[0] and area_size[1] < reduced.size[1]:
+            last_pixel_box = (area_box[2], area_box[3], im.size[0], im.size[1])
+            last_pixel = im.resize((1, 1), Image.BOX, last_pixel_box)
+            reference.paste(last_pixel, area_size)
+
+        self.assert_compare_images(reduced, reference, average_diff, max_diff)
+
+    def assert_compare_images(self, a, b, max_average_diff, max_diff=255):
+        self.assertEqual(a.mode, b.mode, "got mode %r, expected %r" % (a.mode, b.mode))
+        self.assertEqual(a.size, b.size, "got size %r, expected %r" % (a.size, b.size))
+
+        a, b = convert_to_comparable(a, b)
+
+        bands = ImageMode.getmode(a.mode).bands
+        for band, ach, bch in zip(bands, a.split(), b.split()):
+            ch_diff = ImageMath.eval("convert(abs(a - b), 'L')", a=ach, b=bch)
+            ch_hist = ch_diff.histogram()
+
+            average_diff = sum(i * num for i, num in enumerate(ch_hist)) / float(
+                a.size[0] * a.size[1]
+            )
+            self.assertGreaterEqual(
+                max_average_diff,
+                average_diff,
+                (
+                    "average pixel value difference {:.4f} > expected {:.4f} "
+                    "for '{}' band"
+                ).format(average_diff, max_average_diff, band),
+            )
+
+            last_diff = [i for i, num in enumerate(ch_hist) if num > 0][-1]
+            self.assertGreaterEqual(
+                max_diff,
+                last_diff,
+                "max pixel value difference {} > expected {} for '{}' band".format(
+                    last_diff, max_diff, band
+                ),
+            )
+
+    def test_mode_L(self):
+        im = self.get_image("L")
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor)
+            self.compare_reduce_with_box(im, factor)
+
+    def test_mode_LA(self):
+        im = self.get_image("LA")
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor, 0.8, 5)
+
+        # With opaque alpha, an error should be way smaller.
+        im.putalpha(Image.new("L", im.size, 255))
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor)
+            self.compare_reduce_with_box(im, factor)
+
+    def test_mode_La(self):
+        im = self.get_image("La")
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor)
+            self.compare_reduce_with_box(im, factor)
+
+    def test_mode_RGB(self):
+        im = self.get_image("RGB")
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor)
+            self.compare_reduce_with_box(im, factor)
+
+    def test_mode_RGBA(self):
+        im = self.get_image("RGBA")
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor, 0.8, 5)
+
+        # With opaque alpha, an error should be way smaller.
+        im.putalpha(Image.new("L", im.size, 255))
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor)
+            self.compare_reduce_with_box(im, factor)
+
+    def test_mode_RGBa(self):
+        im = self.get_image("RGBa")
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor)
+            self.compare_reduce_with_box(im, factor)
+
+    def test_mode_I(self):
+        im = self.get_image("I")
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor)
+            self.compare_reduce_with_box(im, factor)
+
+    def test_mode_F(self):
+        im = self.get_image("F")
+        for factor in self.remarkable_factors:
+            self.compare_reduce_with_reference(im, factor, 0, 0)
+            self.compare_reduce_with_box(im, factor)
diff --git a/setup.py b/setup.py
index f5024990f..fe024ceae 100755
--- a/setup.py
+++ b/setup.py
@@ -54,6 +54,7 @@ _LIB_IMAGING = (
     "Access",
     "AlphaComposite",
     "Resample",
+    "Reduce",
     "Bands",
     "BcnDecode",
     "BitDecode",
diff --git a/src/PIL/Image.py b/src/PIL/Image.py
index a636a2bbe..5a53b2ad6 100644
--- a/src/PIL/Image.py
+++ b/src/PIL/Image.py
@@ -1777,9 +1777,9 @@ class Image:
            If the image has mode "1" or "P", it is
            always set to :py:attr:`PIL.Image.NEAREST`.
            See: :ref:`concept-filters`.
-        :param box: An optional 4-tuple of floats giving the region
-           of the source image which should be scaled.
-           The values should be within (0, 0, width, height) rectangle.
+        :param box: An optional 4-tuple of floats providing
+           the source image region to be scaled.
+           The values must be within (0, 0, width, height) rectangle.
            If omitted or None, the entire source is used.
         :returns: An :py:class:`~PIL.Image.Image` object.
         """
@@ -1824,6 +1824,39 @@ class Image:
 
         return self._new(self.im.resize(size, resample, box))
 
+    def reduce(self, factor, box=None):
+        """
+        Returns a copy of the image reduced by `factor` times.
+        If the size of the image is not dividable by the `factor`,
+        the resulting size will be rounded up.
+
+        :param factor: A greater than 0 integer or tuple of two integers
+           for width and height separately.
+        :param box: An optional 4-tuple of ints providing
+           the source image region to be reduced.
+           The values must be within (0, 0, width, height) rectangle.
+           If omitted or None, the entire source is used.
+        """
+        if not isinstance(factor, (list, tuple)):
+            factor = (factor, factor)
+
+        if box is None:
+            box = (0, 0) + self.size
+        else:
+            box = tuple(box)
+
+        if factor == (1, 1) and box == (0, 0) + self.size:
+            return self.copy()
+
+        if self.mode in ["LA", "RGBA"]:
+            im = self.convert(self.mode[:-1] + "a")
+            im = im.reduce(factor, box)
+            return im.convert(self.mode)
+
+        self.load()
+
+        return self._new(self.im.reduce(factor, box))
+
     def rotate(
         self,
         angle,
diff --git a/src/_imaging.c b/src/_imaging.c
index b0cbfc665..a55ebb64b 100644
--- a/src/_imaging.c
+++ b/src/_imaging.c
@@ -1787,6 +1787,51 @@ _resize(ImagingObject* self, PyObject* args)
     return PyImagingNew(imOut);
 }
 
+static PyObject*
+_reduce(ImagingObject* self, PyObject* args)
+{
+    Imaging imIn;
+    Imaging imOut;
+
+    int xscale, yscale;
+    int box[4] = {0, 0, 0, 0};
+
+    imIn = self->image;
+    box[2] = imIn->xsize;
+    box[3] = imIn->ysize;
+
+    if (!PyArg_ParseTuple(args, "(ii)|(iiii)", &xscale, &yscale,
+                          &box[0], &box[1], &box[2], &box[3]))
+        return NULL;
+
+    if (xscale < 1 || yscale < 1) {
+        return ImagingError_ValueError("scale must be > 0");
+    }
+
+    if (box[0] < 0 || box[1] < 0) {
+        return ImagingError_ValueError("box offset can't be negative");
+    }
+
+    if (box[2] > imIn->xsize || box[3] > imIn->ysize) {
+        return ImagingError_ValueError("box can't exceed original image size");
+    }
+
+    if (box[2] <= box[0] || box[3] <= box[1]) {
+        return ImagingError_ValueError("box can't be empty");
+    }
+
+    if (xscale == 1 && yscale == 1) {
+        imOut = ImagingCrop(imIn, box[0], box[1], box[2], box[3]);
+    } else {
+        // Change box format: (left, top, width, height)
+        box[2] -= box[0];
+        box[3] -= box[1];
+        imOut = ImagingReduce(imIn, xscale, yscale, box);
+    }
+
+    return PyImagingNew(imOut);
+}
+
 
 #define IS_RGB(mode)\
     (!strcmp(mode, "RGB") || !strcmp(mode, "RGBA") || !strcmp(mode, "RGBX"))
@@ -1801,7 +1846,7 @@ im_setmode(ImagingObject* self, PyObject* args)
     char* mode;
     Py_ssize_t modelen;
     if (!PyArg_ParseTuple(args, "s#:setmode", &mode, &modelen))
-    return NULL;
+        return NULL;
 
     im = self->image;
 
@@ -3235,6 +3280,7 @@ static struct PyMethodDef methods[] = {
     {"rankfilter", (PyCFunction)_rankfilter, 1},
 #endif
     {"resize", (PyCFunction)_resize, 1},
+    {"reduce", (PyCFunction)_reduce, 1},
     {"transpose", (PyCFunction)_transpose, 1},
     {"transform2", (PyCFunction)_transform2, 1},
 
diff --git a/src/libImaging/Imaging.h b/src/libImaging/Imaging.h
index 25c15e758..71dc9c003 100644
--- a/src/libImaging/Imaging.h
+++ b/src/libImaging/Imaging.h
@@ -313,6 +313,7 @@ extern Imaging ImagingRotate270(Imaging imOut, Imaging imIn);
 extern Imaging ImagingTranspose(Imaging imOut, Imaging imIn);
 extern Imaging ImagingTransverse(Imaging imOut, Imaging imIn);
 extern Imaging ImagingResample(Imaging imIn, int xsize, int ysize, int filter, float box[4]);
+extern Imaging ImagingReduce(Imaging imIn, int xscale, int yscale, int box[4]);
 extern Imaging ImagingTransform(
     Imaging imOut, Imaging imIn, int method, int x0, int y0, int x1, int y1,
     double *a, int filter, int fill);
diff --git a/src/libImaging/Reduce.c b/src/libImaging/Reduce.c
new file mode 100644
index 000000000..d6ef92f5b
--- /dev/null
+++ b/src/libImaging/Reduce.c
@@ -0,0 +1,1438 @@
+#include "Imaging.h"
+
+#include <math.h>
+
+#define ROUND_UP(f) ((int) ((f) >= 0.0 ? (f) + 0.5F : (f) - 0.5F))
+
+
+UINT32
+division_UINT32(int divider, int result_bits)
+{
+    UINT32 max_dividend = (1 << result_bits) * divider;
+    float max_int = (1 << 30) * 4.0;
+    return (UINT32) (max_int / max_dividend);
+}
+
+
+void
+ImagingReduceNxN(Imaging imOut, Imaging imIn, int box[4], int xscale, int yscale)
+{
+    /* The most general implementation for any xscale and yscale
+    */
+    int x, y, xx, yy;
+    UINT32 multiplier = division_UINT32(yscale * xscale, 8);
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy_from = box[1] + y*yscale;
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx_from = box[0] + x*xscale;
+                UINT32 ss = amend;
+                for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                    UINT8 *line0 = (UINT8 *)imIn->image8[yy];
+                    UINT8 *line1 = (UINT8 *)imIn->image8[yy + 1];
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss += line0[xx + 0] + line0[xx + 1] +
+                              line1[xx + 0] + line1[xx + 1];
+                    }
+                    if (xscale & 0x01) {
+                        ss += line0[xx + 0] + line1[xx + 0];
+                    }
+                }
+                if (yscale & 0x01) {
+                    UINT8 *line = (UINT8 *)imIn->image8[yy];
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss += line[xx + 0] + line[xx + 1];
+                    }
+                    if (xscale & 0x01) {
+                        ss += line[xx + 0];
+                    }
+                }
+                imOut->image8[y][x] = (ss * multiplier) >> 24;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy_from =  box[1] + y*yscale;
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx_from = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss3 = amend;
+                    for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                        UINT8 *line0 = (UINT8 *)imIn->image[yy];
+                        UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+                        for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                            ss0 += line0[xx*4 + 0] + line0[xx*4 + 4] +
+                                   line1[xx*4 + 0] + line1[xx*4 + 4];
+                            ss3 += line0[xx*4 + 3] + line0[xx*4 + 7] +
+                                   line1[xx*4 + 3] + line1[xx*4 + 7];
+                        }
+                        if (xscale & 0x01) {
+                            ss0 += line0[xx*4 + 0] + line1[xx*4 + 0];
+                            ss3 += line0[xx*4 + 3] + line1[xx*4 + 3];
+                        }
+                    }
+                    if (yscale & 0x01) {
+                        UINT8 *line = (UINT8 *)imIn->image[yy];
+                        for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                            ss0 += line[xx*4 + 0] + line[xx*4 + 4];
+                            ss3 += line[xx*4 + 3] + line[xx*4 + 7];
+                        }
+                        if (xscale & 0x01) {
+                            ss0 += line[xx*4 + 0];
+                            ss3 += line[xx*4 + 3];
+                        }
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, 0,
+                        0, (ss3 * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx_from = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss1 = amend, ss2 = amend;
+                    for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                        UINT8 *line0 = (UINT8 *)imIn->image[yy];
+                        UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+                        for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                            ss0 += line0[xx*4 + 0] + line0[xx*4 + 4] +
+                                   line1[xx*4 + 0] + line1[xx*4 + 4];
+                            ss1 += line0[xx*4 + 1] + line0[xx*4 + 5] +
+                                   line1[xx*4 + 1] + line1[xx*4 + 5];
+                            ss2 += line0[xx*4 + 2] + line0[xx*4 + 6] +
+                                   line1[xx*4 + 2] + line1[xx*4 + 6];
+                        }
+                        if (xscale & 0x01) {
+                            ss0 += line0[xx*4 + 0] + line1[xx*4 + 0];
+                            ss1 += line0[xx*4 + 1] + line1[xx*4 + 1];
+                            ss2 += line0[xx*4 + 2] + line1[xx*4 + 2];
+                        }
+                    }
+                    if (yscale & 0x01) {
+                        UINT8 *line = (UINT8 *)imIn->image[yy];
+                        for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                            ss0 += line[xx*4 + 0] + line[xx*4 + 4];
+                            ss1 += line[xx*4 + 1] + line[xx*4 + 5];
+                            ss2 += line[xx*4 + 2] + line[xx*4 + 6];
+                        }
+                        if (xscale & 0x01) {
+                            ss0 += line[xx*4 + 0];
+                            ss1 += line[xx*4 + 1];
+                            ss2 += line[xx*4 + 2];
+                        }
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                        (ss2 * multiplier) >> 24, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx_from = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss1 = amend, ss2 = amend, ss3 = amend;
+                    for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                        UINT8 *line0 = (UINT8 *)imIn->image[yy];
+                        UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+                        for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                            ss0 += line0[xx*4 + 0] + line0[xx*4 + 4] +
+                                   line1[xx*4 + 0] + line1[xx*4 + 4];
+                            ss1 += line0[xx*4 + 1] + line0[xx*4 + 5] +
+                                   line1[xx*4 + 1] + line1[xx*4 + 5];
+                            ss2 += line0[xx*4 + 2] + line0[xx*4 + 6] +
+                                   line1[xx*4 + 2] + line1[xx*4 + 6];
+                            ss3 += line0[xx*4 + 3] + line0[xx*4 + 7] +
+                                   line1[xx*4 + 3] + line1[xx*4 + 7];
+                        }
+                        if (xscale & 0x01) {
+                            ss0 += line0[xx*4 + 0] + line1[xx*4 + 0];
+                            ss1 += line0[xx*4 + 1] + line1[xx*4 + 1];
+                            ss2 += line0[xx*4 + 2] + line1[xx*4 + 2];
+                            ss3 += line0[xx*4 + 3] + line1[xx*4 + 3];
+                        }
+                    }
+                    if (yscale & 0x01) {
+                        UINT8 *line = (UINT8 *)imIn->image[yy];
+                        for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                            ss0 += line[xx*4 + 0] + line[xx*4 + 4];
+                            ss1 += line[xx*4 + 1] + line[xx*4 + 5];
+                            ss2 += line[xx*4 + 2] + line[xx*4 + 6];
+                            ss3 += line[xx*4 + 3] + line[xx*4 + 7];
+                        }
+                        if (xscale & 0x01) {
+                            ss0 += line[xx*4 + 0];
+                            ss1 += line[xx*4 + 1];
+                            ss2 += line[xx*4 + 2];
+                            ss3 += line[xx*4 + 3];
+                        }
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                        (ss2 * multiplier) >> 24, (ss3 * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduce1xN(Imaging imOut, Imaging imIn, int box[4], int yscale)
+{
+    /* Optimized implementation for xscale = 1.
+    */
+    int x, y, yy;
+    int xscale = 1;
+    UINT32 multiplier = division_UINT32(yscale * xscale, 8);
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy_from = box[1] + y*yscale;
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                UINT32 ss = amend;
+                for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                    UINT8 *line0 = (UINT8 *)imIn->image8[yy];
+                    UINT8 *line1 = (UINT8 *)imIn->image8[yy + 1];
+                    ss += line0[xx + 0] + line1[xx + 0];
+                }
+                if (yscale & 0x01) {
+                    UINT8 *line = (UINT8 *)imIn->image8[yy];
+                    ss += line[xx + 0];
+                }
+                imOut->image8[y][x] = (ss * multiplier) >> 24;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy_from =  box[1] + y*yscale;
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss3 = amend;
+                    for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                        UINT8 *line0 = (UINT8 *)imIn->image[yy];
+                        UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+                        ss0 += line0[xx*4 + 0] + line1[xx*4 + 0];
+                        ss3 += line0[xx*4 + 3] + line1[xx*4 + 3];
+                    }
+                    if (yscale & 0x01) {
+                        UINT8 *line = (UINT8 *)imIn->image[yy];
+                        ss0 += line[xx*4 + 0];
+                        ss3 += line[xx*4 + 3];
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, 0,
+                        0, (ss3 * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss1 = amend, ss2 = amend;
+                    for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                        UINT8 *line0 = (UINT8 *)imIn->image[yy];
+                        UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+                        ss0 += line0[xx*4 + 0] + line1[xx*4 + 0];
+                        ss1 += line0[xx*4 + 1] + line1[xx*4 + 1];
+                        ss2 += line0[xx*4 + 2] + line1[xx*4 + 2];
+                    }
+                    if (yscale & 0x01) {
+                        UINT8 *line = (UINT8 *)imIn->image[yy];
+                        ss0 += line[xx*4 + 0];
+                        ss1 += line[xx*4 + 1];
+                        ss2 += line[xx*4 + 2];
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                        (ss2 * multiplier) >> 24, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss1 = amend, ss2 = amend, ss3 = amend;
+                    for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                        UINT8 *line0 = (UINT8 *)imIn->image[yy];
+                        UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+                        ss0 += line0[xx*4 + 0] + line1[xx*4 + 0];
+                        ss1 += line0[xx*4 + 1] + line1[xx*4 + 1];
+                        ss2 += line0[xx*4 + 2] + line1[xx*4 + 2];
+                        ss3 += line0[xx*4 + 3] + line1[xx*4 + 3];
+                    }
+                    if (yscale & 0x01) {
+                        UINT8 *line = (UINT8 *)imIn->image[yy];
+                        ss0 += line[xx*4 + 0];
+                        ss1 += line[xx*4 + 1];
+                        ss2 += line[xx*4 + 2];
+                        ss3 += line[xx*4 + 3];
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                        (ss2 * multiplier) >> 24, (ss3 * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduceNx1(Imaging imOut, Imaging imIn, int box[4], int xscale)
+{
+    /* Optimized implementation for yscale = 1.
+    */
+    int x, y, xx;
+    int yscale = 1;
+    UINT32 multiplier = division_UINT32(yscale * xscale, 8);
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy =  box[1] + y*yscale;
+            UINT8 *line = (UINT8 *)imIn->image8[yy];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx_from = box[0] + x*xscale;
+                UINT32 ss = amend;
+                for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                    ss += line[xx + 0] + line[xx + 1];
+                }
+                if (xscale & 0x01) {
+                    ss += line[xx + 0];
+                }
+                imOut->image8[y][x] = (ss * multiplier) >> 24;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy =  box[1] + y*yscale;
+            UINT8 *line = (UINT8 *)imIn->image[yy];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx_from = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss3 = amend;
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss0 += line[xx*4 + 0] + line[xx*4 + 4];
+                        ss3 += line[xx*4 + 3] + line[xx*4 + 7];
+                    }
+                    if (xscale & 0x01) {
+                        ss0 += line[xx*4 + 0];
+                        ss3 += line[xx*4 + 3];
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, 0,
+                        0, (ss3 * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx_from = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss1 = amend, ss2 = amend;
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss0 += line[xx*4 + 0] + line[xx*4 + 4];
+                        ss1 += line[xx*4 + 1] + line[xx*4 + 5];
+                        ss2 += line[xx*4 + 2] + line[xx*4 + 6];
+                    }
+                    if (xscale & 0x01) {
+                        ss0 += line[xx*4 + 0];
+                        ss1 += line[xx*4 + 1];
+                        ss2 += line[xx*4 + 2];
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                        (ss2 * multiplier) >> 24, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx_from = box[0] + x*xscale;
+                    UINT32 v;
+                    UINT32 ss0 = amend, ss1 = amend, ss2 = amend, ss3 = amend;
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss0 += line[xx*4 + 0] + line[xx*4 + 4];
+                        ss1 += line[xx*4 + 1] + line[xx*4 + 5];
+                        ss2 += line[xx*4 + 2] + line[xx*4 + 6];
+                        ss3 += line[xx*4 + 3] + line[xx*4 + 7];
+                    }
+                    if (xscale & 0x01) {
+                        ss0 += line[xx*4 + 0];
+                        ss1 += line[xx*4 + 1];
+                        ss2 += line[xx*4 + 2];
+                        ss3 += line[xx*4 + 3];
+                    }
+                    v = MAKE_UINT32(
+                        (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                        (ss2 * multiplier) >> 24, (ss3 * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+void
+ImagingReduce1x2(Imaging imOut, Imaging imIn, int box[4])
+{
+    /* Optimized implementation for xscale = 1 and yscale = 2.
+    */
+    int xscale = 1, yscale = 2;
+    int x, y;
+    UINT32 ss0, ss1, ss2, ss3;
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image8[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image8[yy + 1];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                ss0 = line0[xx + 0] +
+                      line1[xx + 0];
+                imOut->image8[y][x] = (ss0 + amend) >> 1;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] +
+                          line1[xx*4 + 0];
+                    ss3 = line0[xx*4 + 3] +
+                          line1[xx*4 + 3];
+                    v = MAKE_UINT32((ss0 + amend) >> 1, 0,
+                                    0, (ss3 + amend) >> 1);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] +
+                          line1[xx*4 + 0];
+                    ss1 = line0[xx*4 + 1] +
+                          line1[xx*4 + 1];
+                    ss2 = line0[xx*4 + 2] +
+                          line1[xx*4 + 2];
+                    v = MAKE_UINT32((ss0 + amend) >> 1, (ss1 + amend) >> 1,
+                                    (ss2 + amend) >> 1, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] +
+                          line1[xx*4 + 0];
+                    ss1 = line0[xx*4 + 1] +
+                          line1[xx*4 + 1];
+                    ss2 = line0[xx*4 + 2] +
+                          line1[xx*4 + 2];
+                    ss3 = line0[xx*4 + 3] +
+                          line1[xx*4 + 3];
+                    v = MAKE_UINT32((ss0 + amend) >> 1, (ss1 + amend) >> 1,
+                                    (ss2 + amend) >> 1, (ss3 + amend) >> 1);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduce2x1(Imaging imOut, Imaging imIn, int box[4])
+{
+    /* Optimized implementation for xscale = 2 and yscale = 1.
+    */
+    int xscale = 2, yscale = 1;
+    int x, y;
+    UINT32 ss0, ss1, ss2, ss3;
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image8[yy + 0];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                ss0 = line0[xx + 0] + line0[xx + 1];
+                imOut->image8[y][x] = (ss0 + amend) >> 1;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image[yy + 0];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7];
+                    v = MAKE_UINT32((ss0 + amend) >> 1, 0,
+                                    0, (ss3 + amend) >> 1);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6];
+                    v = MAKE_UINT32((ss0 + amend) >> 1, (ss1 + amend) >> 1,
+                                    (ss2 + amend) >> 1, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7];
+                    v = MAKE_UINT32((ss0 + amend) >> 1, (ss1 + amend) >> 1,
+                                    (ss2 + amend) >> 1, (ss3 + amend) >> 1);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduce2x2(Imaging imOut, Imaging imIn, int box[4])
+{
+    /* Optimized implementation for xscale = 2 and yscale = 2.
+    */
+    int xscale = 2, yscale = 2;
+    int x, y;
+    UINT32 ss0, ss1, ss2, ss3;
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image8[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image8[yy + 1];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                ss0 = line0[xx + 0] + line0[xx + 1] +
+                      line1[xx + 0] + line1[xx + 1];
+                imOut->image8[y][x] = (ss0 + amend) >> 2;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] +
+                          line1[xx*4 + 3] + line1[xx*4 + 7];
+                    v = MAKE_UINT32((ss0 + amend) >> 2, 0,
+                                    0, (ss3 + amend) >> 2);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] +
+                          line1[xx*4 + 1] + line1[xx*4 + 5];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] +
+                          line1[xx*4 + 2] + line1[xx*4 + 6];
+                    v = MAKE_UINT32((ss0 + amend) >> 2, (ss1 + amend) >> 2,
+                                    (ss2 + amend) >> 2, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] +
+                          line1[xx*4 + 1] + line1[xx*4 + 5];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] +
+                          line1[xx*4 + 2] + line1[xx*4 + 6];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] +
+                          line1[xx*4 + 3] + line1[xx*4 + 7];
+                    v = MAKE_UINT32((ss0 + amend) >> 2, (ss1 + amend) >> 2,
+                                    (ss2 + amend) >> 2, (ss3 + amend) >> 2);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduce1x3(Imaging imOut, Imaging imIn, int box[4])
+{
+    /* Optimized implementation for xscale = 1 and yscale = 3.
+    */
+    int xscale = 1, yscale = 3;
+    int x, y;
+    UINT32 ss0, ss1, ss2, ss3;
+    UINT32 multiplier = division_UINT32(yscale * xscale, 8);
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image8[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image8[yy + 1];
+            UINT8 *line2 = (UINT8 *)imIn->image8[yy + 2];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                ss0 = line0[xx + 0] +
+                      line1[xx + 0] +
+                      line2[xx + 0];
+                imOut->image8[y][x] = ((ss0 + amend) * multiplier) >> 24;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+            UINT8 *line2 = (UINT8 *)imIn->image[yy + 2];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] +
+                          line1[xx*4 + 0] +
+                          line2[xx*4 + 0];
+                    ss3 = line0[xx*4 + 3] +
+                          line1[xx*4 + 3] +
+                          line2[xx*4 + 3];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, 0,
+                        0, ((ss3 + amend) * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] +
+                          line1[xx*4 + 0] +
+                          line2[xx*4 + 0];
+                    ss1 = line0[xx*4 + 1] +
+                          line1[xx*4 + 1] +
+                          line2[xx*4 + 1];
+                    ss2 = line0[xx*4 + 2] +
+                          line1[xx*4 + 2] +
+                          line2[xx*4 + 2];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, ((ss1 + amend) * multiplier) >> 24,
+                        ((ss2 + amend) * multiplier) >> 24, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] +
+                          line1[xx*4 + 0] +
+                          line2[xx*4 + 0];
+                    ss1 = line0[xx*4 + 1] +
+                          line1[xx*4 + 1] +
+                          line2[xx*4 + 1];
+                    ss2 = line0[xx*4 + 2] +
+                          line1[xx*4 + 2] +
+                          line2[xx*4 + 2];
+                    ss3 = line0[xx*4 + 3] +
+                          line1[xx*4 + 3] +
+                          line2[xx*4 + 3];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, ((ss1 + amend) * multiplier) >> 24,
+                        ((ss2 + amend) * multiplier) >> 24, ((ss3 + amend) * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduce3x1(Imaging imOut, Imaging imIn, int box[4])
+{
+    /* Optimized implementation for xscale = 3 and yscale = 1.
+    */
+    int xscale = 3, yscale = 1;
+    int x, y;
+    UINT32 ss0, ss1, ss2, ss3;
+    UINT32 multiplier = division_UINT32(yscale * xscale, 8);
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image8[yy + 0];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                ss0 = line0[xx + 0] + line0[xx + 1] + line0[xx + 2];
+                imOut->image8[y][x] = ((ss0 + amend) * multiplier) >> 24;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image[yy + 0];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] + line0[xx*4 + 11];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, 0,
+                        0, ((ss3 + amend) * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] + line0[xx*4 + 9];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] + line0[xx*4 + 10];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, ((ss1 + amend) * multiplier) >> 24,
+                        ((ss2 + amend) * multiplier) >> 24, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] + line0[xx*4 + 9];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] + line0[xx*4 + 10];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] + line0[xx*4 + 11];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, ((ss1 + amend) * multiplier) >> 24,
+                        ((ss2 + amend) * multiplier) >> 24, ((ss3 + amend) * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduce3x3(Imaging imOut, Imaging imIn, int box[4])
+{
+    /* Optimized implementation for xscale = 3 and yscale = 3.
+    */
+    int xscale = 3, yscale = 3;
+    int x, y;
+    UINT32 ss0, ss1, ss2, ss3;
+    UINT32 multiplier = division_UINT32(yscale * xscale, 8);
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image8[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image8[yy + 1];
+            UINT8 *line2 = (UINT8 *)imIn->image8[yy + 2];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                ss0 = line0[xx + 0] + line0[xx + 1] + line0[xx + 2] +
+                      line1[xx + 0] + line1[xx + 1] + line1[xx + 2] +
+                      line2[xx + 0] + line2[xx + 1] + line2[xx + 2];
+                imOut->image8[y][x] = ((ss0 + amend) * multiplier) >> 24;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+            UINT8 *line2 = (UINT8 *)imIn->image[yy + 2];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] + line0[xx*4 + 11] +
+                          line1[xx*4 + 3] + line1[xx*4 + 7] + line1[xx*4 + 11] +
+                          line2[xx*4 + 3] + line2[xx*4 + 7] + line2[xx*4 + 11];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, 0,
+                        0, ((ss3 + amend) * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] + line0[xx*4 + 9] +
+                          line1[xx*4 + 1] + line1[xx*4 + 5] + line1[xx*4 + 9] +
+                          line2[xx*4 + 1] + line2[xx*4 + 5] + line2[xx*4 + 9];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] + line0[xx*4 + 10] +
+                          line1[xx*4 + 2] + line1[xx*4 + 6] + line1[xx*4 + 10] +
+                          line2[xx*4 + 2] + line2[xx*4 + 6] + line2[xx*4 + 10];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, ((ss1 + amend) * multiplier) >> 24,
+                        ((ss2 + amend) * multiplier) >> 24, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] + line0[xx*4 + 9] +
+                          line1[xx*4 + 1] + line1[xx*4 + 5] + line1[xx*4 + 9] +
+                          line2[xx*4 + 1] + line2[xx*4 + 5] + line2[xx*4 + 9];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] + line0[xx*4 + 10] +
+                          line1[xx*4 + 2] + line1[xx*4 + 6] + line1[xx*4 + 10] +
+                          line2[xx*4 + 2] + line2[xx*4 + 6] + line2[xx*4 + 10];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] + line0[xx*4 + 11] +
+                          line1[xx*4 + 3] + line1[xx*4 + 7] + line1[xx*4 + 11] +
+                          line2[xx*4 + 3] + line2[xx*4 + 7] + line2[xx*4 + 11];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, ((ss1 + amend) * multiplier) >> 24,
+                        ((ss2 + amend) * multiplier) >> 24, ((ss3 + amend) * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+void
+ImagingReduce4x4(Imaging imOut, Imaging imIn, int box[4])
+{
+    /* Optimized implementation for xscale = 4 and yscale = 4.
+    */
+    int xscale = 4, yscale = 4;
+    int x, y;
+    UINT32 ss0, ss1, ss2, ss3;
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image8[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image8[yy + 1];
+            UINT8 *line2 = (UINT8 *)imIn->image8[yy + 2];
+            UINT8 *line3 = (UINT8 *)imIn->image8[yy + 3];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                ss0 = line0[xx + 0] + line0[xx + 1] + line0[xx + 2] + line0[xx + 3] +
+                      line1[xx + 0] + line1[xx + 1] + line1[xx + 2] + line1[xx + 3] +
+                      line2[xx + 0] + line2[xx + 1] + line2[xx + 2] + line2[xx + 3] +
+                      line3[xx + 0] + line3[xx + 1] + line3[xx + 2] + line3[xx + 3];
+                imOut->image8[y][x] = (ss0 + amend) >> 4;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+            UINT8 *line2 = (UINT8 *)imIn->image[yy + 2];
+            UINT8 *line3 = (UINT8 *)imIn->image[yy + 3];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] + line0[xx*4 + 12] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] + line1[xx*4 + 12] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8] + line2[xx*4 + 12] +
+                          line3[xx*4 + 0] + line3[xx*4 + 4] + line3[xx*4 + 8] + line3[xx*4 + 12];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] + line0[xx*4 + 11] + line0[xx*4 + 15] +
+                          line1[xx*4 + 3] + line1[xx*4 + 7] + line1[xx*4 + 11] + line1[xx*4 + 15] +
+                          line2[xx*4 + 3] + line2[xx*4 + 7] + line2[xx*4 + 11] + line2[xx*4 + 15] +
+                          line3[xx*4 + 3] + line3[xx*4 + 7] + line3[xx*4 + 11] + line3[xx*4 + 15];
+                    v = MAKE_UINT32((ss0 + amend) >> 4, 0,
+                                    0, (ss3 + amend) >> 4);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] + line0[xx*4 + 12] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] + line1[xx*4 + 12] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8] + line2[xx*4 + 12] +
+                          line3[xx*4 + 0] + line3[xx*4 + 4] + line3[xx*4 + 8] + line3[xx*4 + 12];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] + line0[xx*4 + 9] + line0[xx*4 + 13] +
+                          line1[xx*4 + 1] + line1[xx*4 + 5] + line1[xx*4 + 9] + line1[xx*4 + 13] +
+                          line2[xx*4 + 1] + line2[xx*4 + 5] + line2[xx*4 + 9] + line2[xx*4 + 13] +
+                          line3[xx*4 + 1] + line3[xx*4 + 5] + line3[xx*4 + 9] + line3[xx*4 + 13];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] + line0[xx*4 + 10] + line0[xx*4 + 14] +
+                          line1[xx*4 + 2] + line1[xx*4 + 6] + line1[xx*4 + 10] + line1[xx*4 + 14] +
+                          line2[xx*4 + 2] + line2[xx*4 + 6] + line2[xx*4 + 10] + line2[xx*4 + 14] +
+                          line3[xx*4 + 2] + line3[xx*4 + 6] + line3[xx*4 + 10] + line3[xx*4 + 14];
+                    v = MAKE_UINT32((ss0 + amend) >> 4, (ss1 + amend) >> 4,
+                                    (ss2 + amend) >> 4, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] + line0[xx*4 + 12] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] + line1[xx*4 + 12] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8] + line2[xx*4 + 12] +
+                          line3[xx*4 + 0] + line3[xx*4 + 4] + line3[xx*4 + 8] + line3[xx*4 + 12];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] + line0[xx*4 + 9] + line0[xx*4 + 13] +
+                          line1[xx*4 + 1] + line1[xx*4 + 5] + line1[xx*4 + 9] + line1[xx*4 + 13] +
+                          line2[xx*4 + 1] + line2[xx*4 + 5] + line2[xx*4 + 9] + line2[xx*4 + 13] +
+                          line3[xx*4 + 1] + line3[xx*4 + 5] + line3[xx*4 + 9] + line3[xx*4 + 13];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] + line0[xx*4 + 10] + line0[xx*4 + 14] +
+                          line1[xx*4 + 2] + line1[xx*4 + 6] + line1[xx*4 + 10] + line1[xx*4 + 14] +
+                          line2[xx*4 + 2] + line2[xx*4 + 6] + line2[xx*4 + 10] + line2[xx*4 + 14] +
+                          line3[xx*4 + 2] + line3[xx*4 + 6] + line3[xx*4 + 10] + line3[xx*4 + 14];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] + line0[xx*4 + 11] + line0[xx*4 + 15] +
+                          line1[xx*4 + 3] + line1[xx*4 + 7] + line1[xx*4 + 11] + line1[xx*4 + 15] +
+                          line2[xx*4 + 3] + line2[xx*4 + 7] + line2[xx*4 + 11] + line2[xx*4 + 15] +
+                          line3[xx*4 + 3] + line3[xx*4 + 7] + line3[xx*4 + 11] + line3[xx*4 + 15];
+                    v = MAKE_UINT32((ss0 + amend) >> 4, (ss1 + amend) >> 4,
+                                    (ss2 + amend) >> 4, (ss3 + amend) >> 4);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduce5x5(Imaging imOut, Imaging imIn, int box[4])
+{
+    /* Fast special case for xscale = 5 and yscale = 5.
+    */
+    int xscale = 5, yscale = 5;
+    int x, y;
+    UINT32 ss0, ss1, ss2, ss3;
+    UINT32 multiplier = division_UINT32(yscale * xscale, 8);
+    UINT32 amend = yscale * xscale / 2;
+
+    if (imIn->image8) {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image8[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image8[yy + 1];
+            UINT8 *line2 = (UINT8 *)imIn->image8[yy + 2];
+            UINT8 *line3 = (UINT8 *)imIn->image8[yy + 3];
+            UINT8 *line4 = (UINT8 *)imIn->image8[yy + 4];
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx = box[0] + x*xscale;
+                ss0 = line0[xx + 0] + line0[xx + 1] + line0[xx + 2] + line0[xx + 3] + line0[xx + 4] +
+                      line1[xx + 0] + line1[xx + 1] + line1[xx + 2] + line1[xx + 3] + line1[xx + 4] +
+                      line2[xx + 0] + line2[xx + 1] + line2[xx + 2] + line2[xx + 3] + line2[xx + 4] +
+                      line3[xx + 0] + line3[xx + 1] + line3[xx + 2] + line3[xx + 3] + line3[xx + 4] +
+                      line4[xx + 0] + line4[xx + 1] + line4[xx + 2] + line4[xx + 3] + line4[xx + 4];
+                imOut->image8[y][x] = ((ss0 + amend) * multiplier) >> 24;
+            }
+        }
+    } else {
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy = box[1] + y*yscale;
+            UINT8 *line0 = (UINT8 *)imIn->image[yy + 0];
+            UINT8 *line1 = (UINT8 *)imIn->image[yy + 1];
+            UINT8 *line2 = (UINT8 *)imIn->image[yy + 2];
+            UINT8 *line3 = (UINT8 *)imIn->image[yy + 3];
+            UINT8 *line4 = (UINT8 *)imIn->image[yy + 4];
+            if (imIn->bands == 2) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] + line0[xx*4 + 12] + line0[xx*4 + 16] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] + line1[xx*4 + 12] + line1[xx*4 + 16] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8] + line2[xx*4 + 12] + line2[xx*4 + 16] +
+                          line3[xx*4 + 0] + line3[xx*4 + 4] + line3[xx*4 + 8] + line3[xx*4 + 12] + line3[xx*4 + 16] +
+                          line4[xx*4 + 0] + line4[xx*4 + 4] + line4[xx*4 + 8] + line4[xx*4 + 12] + line4[xx*4 + 16];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] + line0[xx*4 + 11] + line0[xx*4 + 15] + line0[xx*4 + 19] +
+                          line1[xx*4 + 3] + line1[xx*4 + 7] + line1[xx*4 + 11] + line1[xx*4 + 15] + line1[xx*4 + 19] +
+                          line2[xx*4 + 3] + line2[xx*4 + 7] + line2[xx*4 + 11] + line2[xx*4 + 15] + line2[xx*4 + 19] +
+                          line3[xx*4 + 3] + line3[xx*4 + 7] + line3[xx*4 + 11] + line3[xx*4 + 15] + line3[xx*4 + 19] +
+                          line4[xx*4 + 3] + line4[xx*4 + 7] + line4[xx*4 + 11] + line4[xx*4 + 15] + line4[xx*4 + 19];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, 0,
+                        0, ((ss3 + amend) * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else if (imIn->bands == 3) {
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] + line0[xx*4 + 12] + line0[xx*4 + 16] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] + line1[xx*4 + 12] + line1[xx*4 + 16] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8] + line2[xx*4 + 12] + line2[xx*4 + 16] +
+                          line3[xx*4 + 0] + line3[xx*4 + 4] + line3[xx*4 + 8] + line3[xx*4 + 12] + line3[xx*4 + 16] +
+                          line4[xx*4 + 0] + line4[xx*4 + 4] + line4[xx*4 + 8] + line4[xx*4 + 12] + line4[xx*4 + 16];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] + line0[xx*4 + 9] + line0[xx*4 + 13] + line0[xx*4 + 17] +
+                          line1[xx*4 + 1] + line1[xx*4 + 5] + line1[xx*4 + 9] + line1[xx*4 + 13] + line1[xx*4 + 17] +
+                          line2[xx*4 + 1] + line2[xx*4 + 5] + line2[xx*4 + 9] + line2[xx*4 + 13] + line2[xx*4 + 17] +
+                          line3[xx*4 + 1] + line3[xx*4 + 5] + line3[xx*4 + 9] + line3[xx*4 + 13] + line3[xx*4 + 17] +
+                          line4[xx*4 + 1] + line4[xx*4 + 5] + line4[xx*4 + 9] + line4[xx*4 + 13] + line4[xx*4 + 17];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] + line0[xx*4 + 10] + line0[xx*4 + 14] + line0[xx*4 + 18] +
+                          line1[xx*4 + 2] + line1[xx*4 + 6] + line1[xx*4 + 10] + line1[xx*4 + 14] + line1[xx*4 + 18] +
+                          line2[xx*4 + 2] + line2[xx*4 + 6] + line2[xx*4 + 10] + line2[xx*4 + 14] + line2[xx*4 + 18] +
+                          line3[xx*4 + 2] + line3[xx*4 + 6] + line3[xx*4 + 10] + line3[xx*4 + 14] + line3[xx*4 + 18] +
+                          line4[xx*4 + 2] + line4[xx*4 + 6] + line4[xx*4 + 10] + line4[xx*4 + 14] + line4[xx*4 + 18];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, ((ss1 + amend) * multiplier) >> 24,
+                        ((ss2 + amend) * multiplier) >> 24, 0);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            } else {  // bands == 4
+                for (x = 0; x < box[2] / xscale; x++) {
+                    int xx = box[0] + x*xscale;
+                    UINT32 v;
+                    ss0 = line0[xx*4 + 0] + line0[xx*4 + 4] + line0[xx*4 + 8] + line0[xx*4 + 12] + line0[xx*4 + 16] +
+                          line1[xx*4 + 0] + line1[xx*4 + 4] + line1[xx*4 + 8] + line1[xx*4 + 12] + line1[xx*4 + 16] +
+                          line2[xx*4 + 0] + line2[xx*4 + 4] + line2[xx*4 + 8] + line2[xx*4 + 12] + line2[xx*4 + 16] +
+                          line3[xx*4 + 0] + line3[xx*4 + 4] + line3[xx*4 + 8] + line3[xx*4 + 12] + line3[xx*4 + 16] +
+                          line4[xx*4 + 0] + line4[xx*4 + 4] + line4[xx*4 + 8] + line4[xx*4 + 12] + line4[xx*4 + 16];
+                    ss1 = line0[xx*4 + 1] + line0[xx*4 + 5] + line0[xx*4 + 9] + line0[xx*4 + 13] + line0[xx*4 + 17] +
+                          line1[xx*4 + 1] + line1[xx*4 + 5] + line1[xx*4 + 9] + line1[xx*4 + 13] + line1[xx*4 + 17] +
+                          line2[xx*4 + 1] + line2[xx*4 + 5] + line2[xx*4 + 9] + line2[xx*4 + 13] + line2[xx*4 + 17] +
+                          line3[xx*4 + 1] + line3[xx*4 + 5] + line3[xx*4 + 9] + line3[xx*4 + 13] + line3[xx*4 + 17] +
+                          line4[xx*4 + 1] + line4[xx*4 + 5] + line4[xx*4 + 9] + line4[xx*4 + 13] + line4[xx*4 + 17];
+                    ss2 = line0[xx*4 + 2] + line0[xx*4 + 6] + line0[xx*4 + 10] + line0[xx*4 + 14] + line0[xx*4 + 18] +
+                          line1[xx*4 + 2] + line1[xx*4 + 6] + line1[xx*4 + 10] + line1[xx*4 + 14] + line1[xx*4 + 18] +
+                          line2[xx*4 + 2] + line2[xx*4 + 6] + line2[xx*4 + 10] + line2[xx*4 + 14] + line2[xx*4 + 18] +
+                          line3[xx*4 + 2] + line3[xx*4 + 6] + line3[xx*4 + 10] + line3[xx*4 + 14] + line3[xx*4 + 18] +
+                          line4[xx*4 + 2] + line4[xx*4 + 6] + line4[xx*4 + 10] + line4[xx*4 + 14] + line4[xx*4 + 18];
+                    ss3 = line0[xx*4 + 3] + line0[xx*4 + 7] + line0[xx*4 + 11] + line0[xx*4 + 15] + line0[xx*4 + 19] +
+                          line1[xx*4 + 3] + line1[xx*4 + 7] + line1[xx*4 + 11] + line1[xx*4 + 15] + line1[xx*4 + 19] +
+                          line2[xx*4 + 3] + line2[xx*4 + 7] + line2[xx*4 + 11] + line2[xx*4 + 15] + line2[xx*4 + 19] +
+                          line3[xx*4 + 3] + line3[xx*4 + 7] + line3[xx*4 + 11] + line3[xx*4 + 15] + line3[xx*4 + 19] +
+                          line4[xx*4 + 3] + line4[xx*4 + 7] + line4[xx*4 + 11] + line4[xx*4 + 15] + line4[xx*4 + 19];
+                    v = MAKE_UINT32(
+                        ((ss0 + amend) * multiplier) >> 24, ((ss1 + amend) * multiplier) >> 24,
+                        ((ss2 + amend) * multiplier) >> 24, ((ss3 + amend) * multiplier) >> 24);
+                    memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+                }
+            }
+        }
+    }
+}
+
+
+void
+ImagingReduceCorners(Imaging imOut, Imaging imIn, int box[4], int xscale, int yscale)
+{
+    /* Fill the last row and the last column for any xscale and yscale.
+    */
+    int x, y, xx, yy;
+
+    if (imIn->image8) {
+        if (box[2] % xscale) {
+            int scale = (box[2] % xscale) * yscale;
+            UINT32 multiplier = division_UINT32(scale, 8);
+            UINT32 amend = scale / 2;
+            for (y = 0; y < box[3] / yscale; y++) {
+                int yy_from = box[1] + y*yscale;
+                UINT32 ss = amend;
+                x = box[2] / xscale;
+
+                for (yy = yy_from; yy < yy_from + yscale; yy++) {
+                    UINT8 *line = (UINT8 *)imIn->image8[yy];
+                    for (xx = box[0] + x*xscale; xx < box[0] + box[2]; xx++) {
+                        ss += line[xx + 0];
+                    }
+                }
+                imOut->image8[y][x] = (ss * multiplier) >> 24;
+            }
+        }
+        if (box[3] % yscale) {
+            int scale = xscale * (box[3] % yscale);
+            UINT32 multiplier = division_UINT32(scale, 8);
+            UINT32 amend = scale / 2;
+            y = box[3] / yscale;
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx_from = box[0] + x*xscale;
+                UINT32 ss = amend;
+                for (yy = box[1] + y*yscale; yy < box[1] + box[3]; yy++) {
+                    UINT8 *line = (UINT8 *)imIn->image8[yy];
+                    for (xx = xx_from; xx < xx_from + xscale; xx++) {
+                        ss += line[xx + 0];
+                    }
+                }
+                imOut->image8[y][x] = (ss * multiplier) >> 24;
+            }
+        }
+        if (box[2] % xscale && box[3] % yscale) {
+            int scale = (box[2] % xscale) * (box[3] % yscale);
+            UINT32 multiplier = division_UINT32(scale, 8);
+            UINT32 amend = scale / 2;
+            UINT32 ss = amend;
+            x = box[2] / xscale;
+            y = box[3] / yscale;
+            for (yy = box[1] + y*yscale; yy < box[1] + box[3]; yy++) {
+                UINT8 *line = (UINT8 *)imIn->image8[yy];
+                for (xx = box[0] + x*xscale; xx < box[0] + box[2]; xx++) {
+                    ss += line[xx + 0];
+                }
+            }
+            imOut->image8[y][x] = (ss * multiplier) >> 24;
+        }
+    } else {
+        if (box[2] % xscale) {
+            int scale = (box[2] % xscale) * yscale;
+            UINT32 multiplier = division_UINT32(scale, 8);
+            UINT32 amend = scale / 2;
+            for (y = 0; y < box[3] / yscale; y++) {
+                int yy_from = box[1] + y*yscale;
+                UINT32 v;
+                UINT32 ss0 = amend, ss1 = amend, ss2 = amend, ss3 = amend;
+                x = box[2] / xscale;
+
+                for (yy = yy_from; yy < yy_from + yscale; yy++) {
+                    UINT8 *line = (UINT8 *)imIn->image[yy];
+                    for (xx = box[0] + x*xscale; xx < box[0] + box[2]; xx++) {
+                        ss0 += line[xx*4 + 0];
+                        ss1 += line[xx*4 + 1];
+                        ss2 += line[xx*4 + 2];
+                        ss3 += line[xx*4 + 3];
+                    }
+                }
+                v = MAKE_UINT32(
+                    (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                    (ss2 * multiplier) >> 24, (ss3 * multiplier) >> 24);
+                memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+            }
+        }
+        if (box[3] % yscale) {
+            int scale = xscale * (box[3] % yscale);
+            UINT32 multiplier = division_UINT32(scale, 8);
+            UINT32 amend = scale / 2;
+            y = box[3] / yscale;
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx_from = box[0] + x*xscale;
+                UINT32 v;
+                UINT32 ss0 = amend, ss1 = amend, ss2 = amend, ss3 = amend;
+                for (yy = box[1] + y*yscale; yy < box[1] + box[3]; yy++) {
+                    UINT8 *line = (UINT8 *)imIn->image[yy];
+                    for (xx = xx_from; xx < xx_from + xscale; xx++) {
+                        ss0 += line[xx*4 + 0];
+                        ss1 += line[xx*4 + 1];
+                        ss2 += line[xx*4 + 2];
+                        ss3 += line[xx*4 + 3];
+                    }
+                }
+                v = MAKE_UINT32(
+                    (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                    (ss2 * multiplier) >> 24, (ss3 * multiplier) >> 24);
+                memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+            }
+        }
+        if (box[2] % xscale && box[3] % yscale) {
+            int scale = (box[2] % xscale) * (box[3] % yscale);
+            UINT32 multiplier = division_UINT32(scale, 8);
+            UINT32 amend = scale / 2;
+            UINT32 v;
+            UINT32 ss0 = amend, ss1 = amend, ss2 = amend, ss3 = amend;
+            x = box[2] / xscale;
+            y = box[3] / yscale;
+            for (yy = box[1] + y*yscale; yy < box[1] + box[3]; yy++) {
+                UINT8 *line = (UINT8 *)imIn->image[yy];
+                for (xx = box[0] + x*xscale; xx < box[0] + box[2]; xx++) {
+                    ss0 += line[xx*4 + 0];
+                    ss1 += line[xx*4 + 1];
+                    ss2 += line[xx*4 + 2];
+                    ss3 += line[xx*4 + 3];
+                }
+            }
+            v = MAKE_UINT32(
+                (ss0 * multiplier) >> 24, (ss1 * multiplier) >> 24,
+                (ss2 * multiplier) >> 24, (ss3 * multiplier) >> 24);
+            memcpy(imOut->image[y] + x * sizeof(v), &v, sizeof(v));
+        }
+    }
+}
+
+
+void
+ImagingReduceNxN_32bpc(Imaging imOut, Imaging imIn, int box[4], int xscale, int yscale)
+{
+    /* The most general implementation for any xscale and yscale
+    */
+    int x, y, xx, yy;
+    double multiplier = 1.0 / (yscale * xscale);
+
+    switch(imIn->type) {
+    case IMAGING_TYPE_INT32:
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy_from = box[1] + y*yscale;
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx_from = box[0] + x*xscale;
+                double ss = 0;
+                for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                    INT32 *line0 = (INT32 *)imIn->image32[yy];
+                    INT32 *line1 = (INT32 *)imIn->image32[yy + 1];
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss += line0[xx + 0] + line0[xx + 1] +
+                              line1[xx + 0] + line1[xx + 1];
+                    }
+                    if (xscale & 0x01) {
+                        ss += line0[xx + 0] + line1[xx + 0];
+                    }
+                }
+                if (yscale & 0x01) {
+                    INT32 *line = (INT32 *)imIn->image32[yy];
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss += line[xx + 0] + line[xx + 1];
+                    }
+                    if (xscale & 0x01) {
+                        ss += line[xx + 0];
+                    }
+                }
+                IMAGING_PIXEL_I(imOut, x, y) = ROUND_UP(ss * multiplier);
+            }
+        }
+        break;
+
+    case IMAGING_TYPE_FLOAT32:
+        for (y = 0; y < box[3] / yscale; y++) {
+            int yy_from = box[1] + y*yscale;
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx_from = box[0] + x*xscale;
+                double ss = 0;
+                for (yy = yy_from; yy < yy_from + yscale - 1; yy += 2) {
+                    FLOAT32 *line0 = (FLOAT32 *)imIn->image32[yy];
+                    FLOAT32 *line1 = (FLOAT32 *)imIn->image32[yy + 1];
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss += line0[xx + 0] + line0[xx + 1] +
+                              line1[xx + 0] + line1[xx + 1];
+                    }
+                    if (xscale & 0x01) {
+                        ss += line0[xx + 0] + line1[xx + 0];
+                    }
+                }
+                if (yscale & 0x01) {
+                    FLOAT32 *line = (FLOAT32 *)imIn->image32[yy];
+                    for (xx = xx_from; xx < xx_from + xscale - 1; xx += 2) {
+                        ss += line[xx + 0] + line[xx + 1];
+                    }
+                    if (xscale & 0x01) {
+                        ss += line[xx + 0];
+                    }
+                }
+                IMAGING_PIXEL_F(imOut, x, y) = ss * multiplier;
+            }
+        }
+        break;
+    }
+}
+
+
+void
+ImagingReduceCorners_32bpc(Imaging imOut, Imaging imIn, int box[4], int xscale, int yscale)
+{
+    /* Fill the last row and the last column for any xscale and yscale.
+    */
+    int x, y, xx, yy;
+
+    switch(imIn->type) {
+    case IMAGING_TYPE_INT32:
+        if (box[2] % xscale) {
+            double multiplier = 1.0 / ((box[2] % xscale) * yscale);
+            for (y = 0; y < box[3] / yscale; y++) {
+                int yy_from = box[1] + y*yscale;
+                double ss = 0;
+                x = box[2] / xscale;
+                for (yy = yy_from; yy < yy_from + yscale; yy++) {
+                    INT32 *line = (INT32 *)imIn->image32[yy];
+                    for (xx = box[0] + x*xscale; xx < box[0] + box[2]; xx++) {
+                        ss += line[xx + 0];
+                    }
+                }
+                IMAGING_PIXEL_I(imOut, x, y) = ROUND_UP(ss * multiplier);
+            }
+        }
+        if (box[3] % yscale) {
+            double multiplier = 1.0 / (xscale * (box[3] % yscale));
+            y = box[3] / yscale;
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx_from = box[0] + x*xscale;
+                double ss = 0;
+                for (yy = box[1] + y*yscale; yy < box[1] + box[3]; yy++) {
+                    INT32 *line = (INT32 *)imIn->image32[yy];
+                    for (xx = xx_from; xx < xx_from + xscale; xx++) {
+                        ss += line[xx + 0];
+                    }
+                }
+                IMAGING_PIXEL_I(imOut, x, y) = ROUND_UP(ss * multiplier);
+            }
+        }
+        if (box[2] % xscale && box[3] % yscale) {
+            double multiplier = 1.0 / ((box[2] % xscale) * (box[3] % yscale));
+            double ss = 0;
+            x = box[2] / xscale;
+            y = box[3] / yscale;
+            for (yy = box[1] + y*yscale; yy < box[1] + box[3]; yy++) {
+                INT32 *line = (INT32 *)imIn->image32[yy];
+                for (xx = box[0] + x*xscale; xx < box[0] + box[2]; xx++) {
+                    ss += line[xx + 0];
+                }
+            }
+            IMAGING_PIXEL_I(imOut, x, y) = ROUND_UP(ss * multiplier);
+        }
+        break;
+
+    case IMAGING_TYPE_FLOAT32:
+        if (box[2] % xscale) {
+            double multiplier = 1.0 / ((box[2] % xscale) * yscale);
+            for (y = 0; y < box[3] / yscale; y++) {
+                int yy_from = box[1] + y*yscale;
+                double ss = 0;
+                x = box[2] / xscale;
+                for (yy = yy_from; yy < yy_from + yscale; yy++) {
+                    FLOAT32 *line = (FLOAT32 *)imIn->image32[yy];
+                    for (xx = box[0] + x*xscale; xx < box[0] + box[2]; xx++) {
+                        ss += line[xx + 0];
+                    }
+                }
+                IMAGING_PIXEL_F(imOut, x, y) = ss * multiplier;
+            }
+        }
+        if (box[3] % yscale) {
+            double multiplier = 1.0 / (xscale * (box[3] % yscale));
+            y = box[3] / yscale;
+            for (x = 0; x < box[2] / xscale; x++) {
+                int xx_from = box[0] + x*xscale;
+                double ss = 0;
+                for (yy = box[1] + y*yscale; yy < box[1] + box[3]; yy++) {
+                    FLOAT32 *line = (FLOAT32 *)imIn->image32[yy];
+                    for (xx = xx_from; xx < xx_from + xscale; xx++) {
+                        ss += line[xx + 0];
+                    }
+                }
+                IMAGING_PIXEL_F(imOut, x, y) = ss * multiplier;
+            }
+        }
+        if (box[2] % xscale && box[3] % yscale) {
+            double multiplier = 1.0 / ((box[2] % xscale) * (box[3] % yscale));
+            double ss = 0;
+            x = box[2] / xscale;
+            y = box[3] / yscale;
+            for (yy = box[1] + y*yscale; yy < box[1] + box[3]; yy++) {
+                FLOAT32 *line = (FLOAT32 *)imIn->image32[yy];
+                for (xx = box[0] + x*xscale; xx < box[0] + box[2]; xx++) {
+                    ss += line[xx + 0];
+                }
+            }
+            IMAGING_PIXEL_F(imOut, x, y) = ss * multiplier;
+        }
+        break;
+    }
+}
+
+
+Imaging
+ImagingReduce(Imaging imIn, int xscale, int yscale, int box[4])
+{
+    ImagingSectionCookie cookie;
+    Imaging imOut = NULL;
+
+    if (strcmp(imIn->mode, "P") == 0 || strcmp(imIn->mode, "1") == 0)
+        return (Imaging) ImagingError_ModeError();
+
+    if (imIn->type == IMAGING_TYPE_SPECIAL)
+        return (Imaging) ImagingError_ModeError();
+
+    imOut = ImagingNewDirty(imIn->mode,
+                            (box[2] + xscale - 1) / xscale,
+                            (box[3] + yscale - 1) / yscale);
+    if ( ! imOut) {
+        return NULL;
+    }
+
+    ImagingSectionEnter(&cookie);
+
+    switch(imIn->type) {
+    case IMAGING_TYPE_UINT8:
+        if (xscale == 1) {
+            if (yscale == 2) {
+                ImagingReduce1x2(imOut, imIn, box);
+            } else if (yscale == 3) {
+                ImagingReduce1x3(imOut, imIn, box);
+            } else {
+                ImagingReduce1xN(imOut, imIn, box, yscale);
+            }
+        } else if (yscale == 1) {
+            if (xscale == 2) {
+                ImagingReduce2x1(imOut, imIn, box);
+            } else if (xscale == 3) {
+                ImagingReduce3x1(imOut, imIn, box);
+            } else {
+                ImagingReduceNx1(imOut, imIn, box, xscale);
+            }
+        } else if (xscale == yscale && xscale <= 5) {
+            if (xscale == 2) {
+                ImagingReduce2x2(imOut, imIn, box);
+            } else if (xscale == 3) {
+                ImagingReduce3x3(imOut, imIn, box);
+            } else if (xscale == 4) {
+                ImagingReduce4x4(imOut, imIn, box);
+            } else {
+                ImagingReduce5x5(imOut, imIn, box);
+            }
+        } else {
+            ImagingReduceNxN(imOut, imIn, box, xscale, yscale);
+        }
+
+        ImagingReduceCorners(imOut, imIn, box, xscale, yscale);
+        break;
+
+    case IMAGING_TYPE_INT32:
+    case IMAGING_TYPE_FLOAT32:
+        ImagingReduceNxN_32bpc(imOut, imIn, box, xscale, yscale);
+
+        ImagingReduceCorners_32bpc(imOut, imIn, box, xscale, yscale);
+        break;
+    }
+
+    ImagingSectionLeave(&cookie);
+
+    return imOut;
+}